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ABSTRACT 


An experimental investigation was conducted to develop a method of predicting 
cylinder indicated torques in a reciprocating engine by measurement of crankshaft angular 
velocity fluctuations. Cylinder indicated pressures were measured for all three cylinders | 
of a two-stroke Diesel engine with pressure transducers. Time-resolved angular position 
was measured at the crankshaft front and at the flywheel. A six degree-of-freedom 
torsional crankshaft model was developed. Two solution methods are described to solve 
the equations of motion: a time-marching ODE solver, and a Finite Element solution in 
the time domain. Using these methods with the measured cylinder torques, the angular 
positions are predicted and compared to measured angular positions for model 
calibration. An inverse solution method was developed to determine the cylinder 
indicated torques from the measured angular position at the crankshaft endpoints. The 
method is theoretically demonstrated to be useful for explicit solutions for two-stroke 
engines up to three cylinders, and four-stroke engines up to four cylinders. Experimental 
results show that the method is useful in predicting cylinder indicated torques from 


angular velocity measurements. 


=> 7 
ll - 
4 , 
4 = o - 
cn) 
ag 
7 
, ss = a a does 





TABLE OF CONTENTS 


Aran TUNG Titec) a) LT aT Os errr 5 eee cc op Se ae'sa deca ve cannes la osleh@UeeetGecaccesccedeneeccaunes | 
Ere 1H aT CLO) Is) irene en ere on cpco coe cosaccenheneecseseetBeetysessdsiesseeeesessenssoe l 

Pe peerless CO) EO TuINNesr ree eo eggs an Seek da sss ooseemer npn Wise acne canacaanneneds 2 
CeO) es N= dU hy less erent eee eee ete acevo asta eastendocnuseesdenvs-cucecvedavesaseeeemmimmtet sso. 5 
rae Ey Nr) I ere re eee year ec ceetae fates stare sdsc-yetecddndsstetters ciiveddceeseaessnantatnns 6 
HA NIG ou ties tt ec TON Aes 1 I aes saan c ccc coca cecuuvencanovebesqsisvenscosacevcssssssescce. 7 
Fen te ls SIG AT snide Verret ern MOR cc. ces 222 < Wace otic sd tn oss ledsbsaveacessaechessdvestverees 7 
Po OAT ONS ey GIN ate ee cee enn Seiavavacacocdcsnedetevececsesedeseoeecsececceatens 8 
CCAR CUIEAIIO NGS av IODEIl Peele TERS i. tiedesecust ctteeancetease. coeeccescisasseuecenees 10 
POs Fe o> ES oN TES NT x TV Dd GD Serer secre cnc... 2. celta secaasnscnscoccettteesss«sssesdacseascssdes 1] 
Ee ENGINE DDE SCRUM MON ete ec sediivesocesaba steers ce dctededeenscceldensec 1] 
J GAN) SOs COIL kN GA) Bl 8 (8k 12) pe i 
rm) (CAT NGO DB arr te satan dca e «480 Es sans ka coi cccsaphsednsdecusevoueedes stance 14 
DD MAGNETIGINDUGBION PROX IMIRIRE RR r-2-.:....1...0s-.se-ossetecacaneccdedetiocsectescesseens 16 
By MODULA TION DOI IN GANA ZER 6 iooeees..cscssyevesennssescstcrteenenadcseennesocececeantes M7 
BES TD V1 Aa earn a bynes saa tu sve ein ceateensoansdconeeteedcentssdecectosacseace 18 
IV. METHODS FOR CYLINDER PRESSURE AND TORQUE PREDICTION .......... 19 
A. PREDICTION OF PHASE DEVIATION FROM CYLINDER TORQUES ......... i 
fe) Tine = aan ima em aaa OU 7. meee co 2cs os. sanedaste< coos cecieeeene-cesetne ache tecesevessveees 19 

Pee Ae we Cis MMMUINC HOME <See es... 2, eter SMM nas sta enec oat. Sanaa ese edaty lds eee 24 

B. PREDICTION OF TORQUES FROM MEASURED SHAFT SPEEDS .............. Zo 
Ke SOLUTIONS len Anne le GU UO emer coee ee othe. 20s dacassncsenae Meet eeteenteee. ©... che: 2S 

2 ANCE TPO labtop @ tee ait eran te mene tere 23 2565.0 0025 2a salnsese Meee dete <6 od ss oo 0 28 

See ore nal mtemime Dy loco wl O lien lene S GON, cates 2--20ce02.--+..ncee cee cose pie. -6 <5 ese Oy 

Ase bh egtes ter ll DNV Re Meany 21. 10) , 9 AURIS ei en nee ee eee Zo 

Re ES ee SO PCa IN TO oe cece cece ce. ceenge ede a sonsincn Anes seacesa¢vesesesvedbelesi teers 30 
eats es SLB ALE Sie oe ee meeteea et aE ie 85  , Rer ee eRe reese ala val aoscien Soe. She cue tn ee etetousibtoviussaesee 33 
pow CALIB RAIONIOE Ne RW NIOD BU ovccncyeccsccastis-sesssiescsssdvevenvossharesedescredene 33 
B. COMPARISON OF PREDICTED AND MEASURED TORQUES ..........0.00000- By 
VI. SUMMARY, CONCLUSIONS, AND RECOMMENDATIONS. ...................::.0cee0ee 4] 
en ) VIVE RO eres ee ee eee MMMM So 2h cas «045 se0s.s-san<cecccessceceedsvenscecercusescessvee 4] 
ea IN CUS 1G sie cere fad 002506 Rs Gua esas ci sedsse~ sec -sseeccdaned caveeeacecciceesceeeee 4] 
Oe REC OMNI Be INMD MN INS eet sooo en c= ceeds ea cea vo bcc vnissntectesegetinss cossteevesveeceeseuss 42 
APPENDIX A] GEOMETRY OF ROTATING COMPONENTS .....c1........0.ccccccsteeecoses 45 
APPENDIX B. GEOMETRY OF RECIPROCATING COMPONENTS ...................04. 55 
APPENDIX C. NATURAL FREQUENCY AND MODAL ANALYSIS .................... oe 
ewratate INLET TAT A) NT Sy ce 2 e.vne one c sees cde sense tctenscosseesdonn ones 63 
Pesteg ite IN NEAL TNA IE iss ODT Nett eee mee ec aces vs cecsdavevoadvsssos.cdaccnaeesseederosssecnveendeacs 83 
Pees eg IN LEY se LIN a Trafic o..c0c2ccceace ence edeccecesscsad¥acedteeedeceeeeesesatense 103 
PU STO) Fe aba Fe Fe Fe dle ea re cee 2 os ce cade «cc acdna Jadu .ovesuiennanoesbis soacosqeBeaeeceosareetenes ctbevers 105 
STILUA cede ee) USA elites NO MDW) IML Sy eens tet teh oy insdseaaaclaveecacascceteesscéenneccedees vesovevessteaoeeess 109 





m 


: 
7 

on 
Oo 


Gare 


OU 


a | > 99 aT = 
7 ee 
a 


— 
‘tiem 


NOMENCLATURE 


Description 
Piston cross-sectional area 


Cylinder bore 

Torsional damping 

Damping matrix 

Diameter 

Crankpin diameter 

Main journal bearing diameter 
Friction factor 

Connecting rod force 

Net cylinder force from indicated pressure 
Reaction force 

Element torque vector 

Finite Element Method torque vector 
Gravitational acceleration 

Finite Element duration 

Average bearing clearance 


Shape functions 

Weighted average of residual 

Identity matrix 

Rotating mass polar moment of inertia 
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Connecting rod rotating inertia 
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Finite Element Method matrix 
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Crankshaft rotational velocity 
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Piston linear position 
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J. INTRODUCTION 


A. MOTIVATION 


Reciprocating engines require a maintenance system to ensure readiness 
throughout operational life. Currently, particularly in military uses, this 1s accomplished 
with a regularly scheduled maintenance (RSM) system, where parts are checked, and 
replaced if necessary, and components are serviced at intervals based on_ historical 
expected failure rates. Therefore, maintenance is normally performed well before it is 
actually necessary, because the true condition of the engine is unknown. Significant 
savings could be realized with the use of a condition based maintenance (CBM) system. 
Such a system requires a means to monitor engine health during operation, or with 
operational tests that do not require significant work on the engine. 

Several classes of faults that occur in Diesel engines can be detected and localized 
by measurement of individual cylinder firing pressures. Examples include loss of 
compression ratio due to cylinder leaks, and improper combustion of fuel due to injector 
problems. Monitoring cylinder firing pressure is an excellent means of condition based 
maintenance. However, due to the harsh environment in the cylinders, the pressure 
transducers required are very expensive and short-lived. While direct measurement of 
cylinder pressures for performance monitoring is feasible and sometimes used in an 
Operational engine, it is expensive. Of course, a possible solution to this problem would 
be the development of cheap, reliable pressure transducers for use in operational engines. 


Barring this, an alternative solution is the use of indirect methods for estimating cylinder 


pressures, such as the measurement of crankshaft angular velocity fluctuations along with 


an appropriate scheme for inferring the pressure waveform. 


B. STATE OF THE ART 


The angular velocity of a reciprocating engine contains small fluctuations due to 
the variations of cylinder pressures. In general, the engine speeds up after a cylinder 
fires, then slows down as the next cylinder is compressing in preparation for combustion. 
The flywheel is intended to reduce the magnitude of these oscillations, but they are still 
present and represent a speed variation of several percent, which is a measurable amount. 
A number of researchers have investigated the possibility of predicting the cylinder 
pressure variation by measurement of these small speed oscillations. 

Freestone and Jenkins [Ref 1] measured crankshaft velocity with a proximeter 
mounted at the flywheel ring gear teeth. They developed a lumped crankshaft model, 
using inertial torque to account for the reciprocating piston masses. This model was used 
to develop a calculation of the total gas torque in the engine cylinders as a function of 
crank angle. Noting abnormally low peaks in the pressure waveform and _ their 
corresponding crank angle localized faults in individual cylinders. 

Mauer and Watts [Ref 2] measured angular velocity at both ends of the crankshaft 
by placing a proximeter at the flywheel ring gear teeth and at a corresponding gear 
mounted on the pulley. The phase difference between the two encoders corresponded to 
an instantaneous measurement of the total crankshaft twist, which was considered 
proportional to crankshaft torque. No mathematical model was used, so detection of 
faults was realized by comparing the measured signal to a signal recorded on a healthy 


engine. As expected, the twist signal had peaks that varied depending on which cylinder 
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was currently firing, because cylinders farther from the flywheel cause a greater twist for 
the same combustion pressure. Mauer continued this work [Ref 3], developing a lumped- 
parameter engine model to isolate cylinder specific torque. In this method, he computed 
the cylinder specific torque, defined as the integration of the total crankshaft torque from 
TDC of the cylinder in question to TDC of the next firing cylinder. In a healthy engine, 
cylinder torques will all be close to the mean, but a faulty cylinder will show a reduced 
specific torque compared to the other cylinders. 

Citron, et al. [Ref 4] used a four degree-of-freedom model of the engine-drivetrain 
system that differentiated between the flywheel and the engine. Reciprocating masses 
were accounted for by an inertial torque component and the crankshaft speed was 
measured with a proximeter at the flywheel ring gear teeth. Total cylinder gas pressures 
were reconstructed by solving the equations of motion. Individual cylinder pressures 
were inferred by assuming that the majority of the net torque at a particular point was due 
to the cylinder undergoing the power stroke. 

Connolly and Yagle [Ref 5] used a lumped engine model, assuming the total 
inertia of the crankshaft components as a single mass. An inertial torque accounted for 
reciprocating masses and the angular velocity was measured with a proximeter at the 
flywheel ring gear. A nonlinear differential equation relating combustion pressure to 
angular velocity, derived from the torque balance equation, was reformulated to a linear 
first-order differential equation relating pressure to the square of the angular velocity. 
Connolly revisited the issue [Ref 6] to reconstruct cyclic pressure variability from the 


crankshaft angular velocity. 


Lim et al. [Ref 7] predicted cylinder pressures in a four-cylinder four-stroke spark 
ignition engine by making several assumptions about the cylinder pressure as a function 
of the crank angle. The intake and exhaust strokes were assumed reference values (intake 
manifold pressure and exhaust backpressure) and the compression stroke was estimated 
as a polytropic process for each cylinder. From knowledge of the crank angle and the 
firing order, power stroke pressures were estimated for each cylinder from the measured 
angular velocity and the known load torque. The method implied a lumped crankshaft 
model. 

lida et al. [Ref 8] measured angular velocity with a proximeter at the flywheel, 
and included a correction for tooth-to-tooth variation, determined by measurement of the 
flywheel at a constant rotational speed. A lumped engine model was used, in which an 
equation related the total engine inertia and rotational acceleration to the composite 
torque applied to the crankshaft. Integration of this equation over a cycle yielded a 
relation to determine Indicated Mean Effective Pressure (IMEP) from the total engine 
inertia and the square of the change in the angular velocity. 

Taraza [Ref 9] developed a linear high-fidelity model of a multicylinder engine. 
Using this model, angular velocity and angular deflection were predicted for the front of 
the crankshaft and compared to measured values for an inline, four-stroke, four-cylinder 
Diesel engine, in order to verify the model parameters. He then determined harmonic 
orders for the crankshaft model, and conducted experimental measurements at these 
speeds. His results show good agreement between the measured and predicted harmonic 
order amplitudes. His conclusion was that the measured amplitudes of certain harmonic 


orders of angular motion could be used to determine engine mean indicated pressure. 


Additionally, his work demonstrated the usefulness of a high-fidelity crankshaft torsional 
model. 

Additional work on this subject [Refs 10-16] generally used a lumped crankshaft 
model with angular velocity measured at one point. 

Previous work at NPS by Bell {Ref 17] and Hudson [Ref 18], on the same engine 
used in this study, demonstrated that there is information present in the crankshaft 
rotational speed of a reciprocating engine. Hudson developed a high-fidelity model of 
the engine crankshaft, then used measured pressures to predict speed fluctuations for 
comparison with actual speed fluctuations at the crankshaft nose. Due to noise from 
vibration in the optical encoder mounting, he was unable to show good agreement for 
speed fluctuations between measured results and the model. 

To the best of the author’s knowledge, no research has been reported using a 
high-fidelity torsional model to determine explicit cylinder-specific torques throughout a 
representative cycle. An engine crankshaft displays torsion that varies during a cycle, 
and along its length. This twist absorbs rotational energy that is later released when the 
twist relaxes. Previous models that consider the crankshaft as one rotating element 
neglect the effect this twist produces on the torques of individual cylinders. But when the 
intent is to localize engine faults, it is imperative that these differences between cylinders 


are considered. 


Cc: OBJECTIVES 


The primary objective of this study is to develop a method of determining 
individual cylinder gas torques from measured time-resolved angular positions at the two 


endpoints of the crankshaft, using a high-fidelity torsional model of the crankshaft. 
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Calibrated parameters will be determined for the torsional model from calculated values 
and experimental data. 

Another objective is to develop numerical solution techniques for solving the 
equations of motion in both directions. Specifically, a method for direct integration of 
the differential equations of motion will be formulated that sets cyclic boundary 
conditions in the time domain. The reverse method, to determine cylinder gas torques 
from tine-resolved angular position data, must use data from only two measurement 


points instead of all degrees of freedom. 


D. ORGANIZATION 


Chapter II describes the development and calibration of the crankshaft torsional 
model. The equations of motion for the crankshaft are derived and presented. 

Chapter III describes the experimental apparatus used to collect data for this 
study. Specifications for the test engine and diagrams for the instrumentation are | 
presented. A summary of the various engine operating conditions used for the study 1s 
included. 

In Chapter IV, an explanation of the numerical methods used to solve the 
equations of motion is described. This will include numerical methods for model 
calibration as well as cylinder gas torque prediction. 

Chapter V shows the results of the engine test runs and data analysis that will 
support the thesis concept. 

Chapter VI summarizes the study, presents conclusions, and _ lists 


recommendations for further research in this area. 


Il. CRANKSHAFT TORSIONAL MODEL 


A. PHYSICAL SYSTEM 

The engine used in this research is a Detroit Diesel (3-53 Series) three-cylinder 
two-stroke Diesel. The crankshaft (Figure 1) 1s supported by four main journal bearings, 
and includes counterweight lobes on four of the six crankwebs. The front of the engine is 
to the left in the diagram. A press fitted gear at the crankshaft nose drives the oi] pump, 


and auxiliary loads are driven by the timing gear, located just forward of the flywheel 


(not shown). 
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Figure 1. Crankshaft. From Ref [19]. 


An idealized mass-elastic torsional model is used to mathematically describe the 
angular motion of the crankshaft (Figure 2). This model, originally developed by Hudson 
[Ref 18], has been refined for the present study. Specifically, additional load torques 
were added to the model to account for the effect of the oil pump and the auxiliary loads, 
constant parasitic force was used to model the piston ring friction, and the effect of 


reciprocating torques was added. The model consists of six concentrated masses 
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connected by massless idealized shafting. The mass concentrations are centered at the 
optical encoder mounting, the three crankpins, the flywheel, and the dynamometer rotor. 
Torsional rigidity and damping are indicated by K and C, respectively. Gas torques and 


load torques are applied at the mass concentrations, as indicated in parenthesis. 
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Figure 2. Crankshaft Torsional Model. After Ref [18]. 


B. EQUATIONS OF MOTION 


Using the model from Figure 2, a set of six second order differential equations are 
developed to describe the rotational dynamics of the crankshaft. Since no part of the 
crankshaft is fixed, the model requires six separate angular position indications. These 
are the crank angles at each of the lumped mass points designated in the model, and they 


are designated 0, through 6. 


JO: + C1281 - @2) + Ki(01 — 82) = —Tpmp 


(1) 


(J2+J,,,,)02+ C1{G2—@) + Ki(@2— 01) + CA G2 - 3) + Ki. - 3) + C22 =T,,, (H+ T,,,(t)-T,, (2) 
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J 606+ Cs6( O6 = @s) + K5( 06 — @s) = —Ttoad (6) 
For simplicity, the six equations can be combined into one matrix equation: 
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G CALCULATION OF MODEL PARAMETERS 


The moments of inertia and the torsional rigidities can be analytically calculated 
for the model, as described in detail in Appendices A and B. It is assumed that the values 


for damping will be very low. Final values for the model are listed in Table 1. 


Table 1. Model Parameter Values 
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Values for the friction and auxiliary loads are more difficult to determine 










analytically. Additionally, they will vary depending on the load and speed of the engine. 
Appendix A contains a theoretical analysis of friction and load torques, and this analysis 
was used to formulate an estimate of the expected magnitudes of friction and load 


torques. The values actually used in the model are listed in Table 2. 


Table 2. Model Friction and Load Values 
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Values for the nonlinear model parameters (Trec, Jrec, and Tey) are determined as 














functions of 8 or t. Their derivation is described in detail in Appendices A and B. 


Il. EXPERIMENTAL METHODS 


A. ENGINE DESCRIPTION 


The Engine used in this research was a Detroit Diesel 3-53 Series engine, with 
characteristics listed in Table 3. For this study the front of the engine is designated as the 
end where the pulley would be located; the rear is the flywheel end. The cylinders are 
numbered consecutively from front to rear, so that cylinder #1 is the farthest from the 
flywheel. Cylinders are naturally aspirated; a roots blower provides a positive crankcase 
pressure that is proportional to engine speed [Ref 18]. The engine is considered to be a 


typical example of a Diesel engine. 


Table 3. Engine Characteristics. From Ref [19] 


Model 5033-5001N 

Type 

Number of Cylinders 

Number of Main Bearings 
4 


Firing Order 3-2; Clockwise Rotation 














Exhaust Valves per Cylinder 


Displacement per Cylinder 


Compression Ratio 
Bore 
Stroke 
Max Rotation Speed 





The engine has been slightly modified. The front-end pulley was removed for 
mounting of the optical encoder to the crankshaft, and the alternator was removed (Figure 
3). The engine was mounted on a Superflow engine test stand, and was loaded by an SF- 


901 Water Brake Dynamometer (Figure 4). 
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Figure 3. Engine Test Stand (Front View) 
Figure 4. Engine Test Stand 
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B. ENGINE CYCLE ANALYZER (ECA) 


The Superflow Engine Cycle Analyzer (Figure 5) is a PC based data acquisition 
system. A sensor interface collects an engine load signal from the dynamometer, a crank 
angle and TDC signal from the optical encoder, and pressure signals from the 
piezoelectric pressure transducers mounted in the glow plug sockets for each cylinder. 
This oneton is passed to a data acquisition computer, which is used to store and 
display the pressure data. Raw pressure data are collected and phase-lock ensemble 


averaged over |1 cycles. then used in the numerical analysis programs. [Ref 20] 
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Figure 5. Instrumentation Schematic 
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C. OPTICAL ENCODER 


A Heidenhain incremental Rotary Encoder [Ref 21] was used to collect time- 
resolved angular position at the crankshaft nose. The encoder consists of a flat optical 
disk, rigidly attached to the rotating shaft, with a specified number of evenly spaced 
windows etched near the perimeter (Figure 6). A signal is generated by photoelectric 
scanning of the disk as it rotates. The output signal is a TTL square wave, where highs 
correspond to a window passing in front of the detector. Measurement of the leading 
edge of the square wave corresponds to a time stamp for a specific angular position. The 
time differences inversely correspond to the average speed of the shaft as it rotates 
through the incremental angle. Encoders with 720 and 3,600 windows were available, 
but in either case 720 counts per revolution were collected, for an angular resolution of 


0.5°. During a run data were collected for 11 cycles at the encoder for a total of 7,920 


time stamps. 





Figure 6. Optical Disk Representation 


The Optical Encoder shaft was mounted to the end of the crankshaft with a 
flexible coupling (Figure 7). The coupling allows for radial and axial vibration of the 
crankshaft that would damage the encoder, because the endplay of the crankshaft exceeds 


the design specifications for the optical encoder without the protective coupling installed. 


Hudson [Ref 18] collected his data with the encoder mounted directly to the crankshaft, 
resulting in extremely noisy data that did not compare favorably with predicted data. 
Additionally, excessive crankshaft radial vibrations damaged several encoders. The 
coupling transmits the angular position of the crankshaft nose to within an accuracy of 
10" (4.85e-05 radians) [Ref 21]. An additional effect of the coupling is a high frequency 
torsional oscillation due to the natural frequency of the coupling/rotor combination. This 
is discussed in detail in Appendix C, and was not a significant problem since the signal of 


interest was at a much lower frequency. 
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Figure 7. Optical Encoder and coupling. From Ref [21] 


The body of the encoder is rigidly mounted to the engine block to minimize noise 
due to vibration (Figure 8). The mounting of the encoder is extremely important. 
Hudson used several different mounting schemes, before settling on the mount used again 


in this study, which works very well to ensure engine vibration does not affect the 


encoder measurements. 
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Figure 8. Optical Encoder Mounting. From Ref [18] 


et OF 
Ste ts Ae — 


D. MAGNETIC INDUCTION PROXIMETER 


A Bentley Nevada 3000 series 190 Proximitor system was used to detect passage 
of the flywheel ring gear teeth. The system consists of a ferromagnetic eddy current 
detector, which outputs a negative voltage that 1s a function of the distance between the 
probe end and a ferromagnetic surface. A TTL conversion circuit triggers a step change 
in voltage when the proximeter output exceeds a certain level, corresponding to a 
Aeon of about '4 inch. The probe was mounted on a bracket fastened across the edge 
of the flywheel access panel, so that it saw the sides of the gear teeth as they passed 
(Figure 9). The output of the circuit was a square wave TTL signal; the leading edge of 
each wave corresponding to the passage of a gear tooth beneath the probe. The TTL 
output signal was sent directly to the MDA for data collection. There are 126 teeth on the 
flywheel ring gear, and during a run data were collected for 63 cycles, for a total of 7938 


time stamps. 
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Figure 9. Proximeter mounting 


E. MODULATION DOMAIN ANALYZER 


A Hewlett Packard 53310A Modulation Domain Analyzer (MDA) was used to 
collect the time stamp data from the optical encoder and the proximeter. In either case, a 
TDC indicator (an output from the ECA) triggered the MDA. It received a TTL signal 
and recorded a time stamp at the leading edge of each wave. The MDA was able to 
collect up to 8,000 data points at a time. A single run collected 11 cycles from the optical 
encoder or 63 cycles from the flywheel proximeter, as described previously. A data 
acquisition computer controlled the MDA and received a transferred file containing the 


time Stamps. 
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F. TEST MATRIX 


A series of data runs were performed to test the validity and consistency of the 
model at varying engine speeds and loads. During a run data were collected for rotational 
speed at the flywheel and crankshaft nose, cylinder pressures for the three cylinders, 
dynamometer load, and atmospheric pressure. A data series was composed of data 
collected during a single run of the engine, at varying loads and speeds, normally taking 
about an hour on a single day. A prefix letter, such as “S,” designated a particular series 
so that comparisons could be restricted to data taken during a single operation of the 
engine. This was intended to eliminate any variations in operation that might take place 
as the engine condition varied over time. Table 4 shows the elements of each data series, 
indicating what variations in load and speed make up each one. Table 5 lists specific 


information for each data run. 


Table 4. Speeds and Loads for Data Series 


DYNAMOMETER LOAD (FT*LBF) 
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Table 5. Data Run Information 


Pressure 
( a Sla) 
Heli — Cal coupling used; slippage of 
encoder shaft 


V 14 Oct 98 14.730 Heidenhain K17 coupling used; TDC lag 
21 Oct 98 14.706 Heidenhain K17 coupling used; no lag 
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IV. METHODS FOR CYLINDER PRESSURE AND TORQUE PREDICTION 

The torsional model may be used to predict the motions of the crankshaft from 
measured applied forces, or it may be used to predict the forces from the measured 
motions. In order to do this, the differential equations of motion must be solved in both 
directions. First, two methods will be described for determining the angular positions 6, 
through 06, given the measured gas torques from the three cylinders. These solution 
methods, referred to as direct integration methods, will be used to test the validity of the 
torsional model and calibrate the parameters. Second, a method will be described for 
determining the individual cylinder gas torques Ti,y through T3,yj, given the time- 
resolved angular position at the two ends of the crankshaft, 6; and @s. This final solution 
method, called the inverse method, will be used for detection of cylinder faults from 


measurements of crankshaft rotational velocity. 


A. PREDICTION OF PHASE DEVIATION FROM CYLINDER TORQUES 


The equations of motion for the crankshaft model constitute a system of non- 
linear, second-order ordinary differential equations (ODEs). Calibration of the values for 
the torsional model will be conducted by first solving the differential equation for {8}, 
given the cylinder indicated torques. The predicted phase deviation and twist determined 
from {8} can be compared to measured values to determine validity of the model. 

I Time-marching O.D.E. method 

The first method is a direct integration of the ODEs in the time domain. This is 
accomplished numerically using a fourth- and fifth-order Runge-Kutte method. First, the 


six second-order ODEs must be converted to 12 first-order ODEs as follows: 
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(8) 


These equations cannot be solved all in one step, however, because they are 
nonlinear due to the dependence of J,-, on 0, and because the T,yjs and Trees are functions 
of time. Instead, a “time marching” method 1s used where the equations are solved over 
small steps and the final condition of each step becomes the initial condition for the 
subsequent step. The values of Tey, Trec, and Jrec can then be approximated as a constant 
value over the step, or as a linear interpolation within the function describing the 
equation. 

This method requires significant computing time. This is because a series of 
“shooting” iterations must be conducted to determine the initial conditions that yield 
cyclic conditions for the representative cycle (1.e., a periodic solution). Since the domain 
is an assumed representative cycle, it follows that the values of 6 and 6 at the end of the 
cycle must match those at the beginning of the cycle. The initial angular velocities 


chosen will have a significant effect on whether or not the solution @ 1s “cyclic.” 
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2 Finite Element method 

Although the integration of the differential equations is an initial value problem, 
because it is assumed to describe a cyclic process the solution of the equation is crown at 
future times. That is, at the end of one representative cycle, we expect all the angular 
positions to be increased by exactly 27 radians. An alternative method for solution 
avoids the shooting iterations by setting “boundary” conditions in time, instead of initial 
conditions. The problem ts then treated as a boundary value problem in time, and a Finite 
Element Method (FEM) is developed to solve a second-order differential equation for 0 
as a function of t. This method is based on Kwon and Bang [Ref 22]. 

The weak formulation of the weighted residual] method 1s used to approximate the 
solution to a second-order matrix differential equation. To accomplish this, the weighted 
average of the residual over the domain is set to zero: 


I= [wf 18 +1C16 +1K}9 -r Jr = {0} (9) 


and then simplified to: 
Uy 
{fle} wl 6 + w[K p Yr = {wr }ar— fol} (10) 

where w is the weighting function, and 0 is a vector corresponding to the angular position 
at the six degrees-of-freedom. 

A Galerkin Finite Element formulation is developed, using the sum of simple 
piecewise linear shape functions to approximate the more complex real function. The 
shape functions used are set up to be | at a node, linearly decreasing to O at the adjacent 


nodes. The value of the function at any point is approximated as a linear combination of 


Z| 
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the values of the two shape functions for the adjacent nodes. Figure 10 shows these 
linear shape functions for a hypothetical element. The function T(t) is approximated as 
FS Hy (t) T(t) + H>(t)T(t,) between the nodes. This becomes the trial function for 
Galerkin’s method, and the test functions are w,;=H)(t) and w);=H>(t). 
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Figure 10. Linear Shape Functions. After Ref [22] 


For a one-dimensional function, the shape functions are determined by: 
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But when used for a 6x6 matrix equation, the test functions must also be expanded into 


the following matrices: 
] | 
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When substituted in Equation (10) and solved over the element from t, to tp, the result is: 


eV =f) (13) 


where 


ere |e pele Hale cl Ha) ac H\\a (14) 


2 
Ty ) 


- 


and 


22 


saa | H h\ 7 
fea ‘ 
Here, the element matrix [K*] is 12x12 and the element torque vector {F°} is 12x1. The 
values of the torque at the two nodes, {T,} and {Tr}, are 6x1 vectors. Equation (13) may 
be solved for the 6x1 vectors {9_} and {9p}, which are the approximate solutions at the 
nodes. The solution used in this study has 720 elements and 721 nodes to solve for one 
representative cycle. The element matrices for each node are assembled into a 
4326x4326 finite element matrix [Ry and the element torque vectors are assembled into 
a 4326x1 finite element torque vector {F°}. Boundary conditions are established by 
defining the value of the 6x1 vector {8} at the end nodes. This is accomplished by 
setting the first and last six lines of the finite element matrix (K*] to identity, and setting 
the first six and last six values in the finite element torque vector {F°} to the boundary 
values. The solution {9} (a 4326x1 vector, the 720 6x1 nodal solutions {0;} stacked 
vertically) is then found from the matrix equation: 
IK" 6 =F" (16) 

This method shows a marked improvement in computing efficiency over the time- 
marching method. In addition to being about three times as fast for each program run, the 
Finite Element Method avoids the shooting iterations required for the time-marching 
method, which had to be repeated as many as five times for each solution. A comparison 
of the Phase Deviation found with both methods is shown in Figure 11. The Phase 


Deviation Plot will be described in detail in the next chapter. 
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Figure 11. Comparison of Solution Methods 
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B. PREDICTION OF TORQUES FROM MEASURED SHAFT SPEEDS 
i Solution of Matrix Equation 

a. Two-Stroke Engines 

Given a properly calibrated model, the time dependent values of gas 
torque for the three cylinders (T]cy], T2cy], and T3cy]) can be determined from the six 
second-order ODEs (Equations 1-6). Figure 12 outlines the soiution method. By 
measurements at the flywheel and crankshaft nose, the values for 0; and 0; can be 
determined for a representative cycle. Numerical differentiation of this data yields the 
velocities and accelerations at these two points. 

With properly calibrated parameters for the torsional model, Equation (6) 
can be solved for Q,., then Equation (5) can be solved for 84, and Equation (1) can be 
solved for 82. Now there are four unknowns left (03, Ticy, Tacy, and T3,y) and three 
equations (Equations 2, 3, and 4). However, because this is a two-stroke engine, each 
cylinder is at reference pressure for one third of each cycle, while the exhaust and intake 
ports are uncovered. Therefore, at any one time during a representative cycle, one of the 
three cylinder torques is known, leaving three unknowns and three equations. For 
example, for the first 120 degrees of the cycle, cylinder #2 ports are open, so the pressure 
in cylinder #2 is reference pressure (See Figure 13). For the first 120 degrees, T3,y, can 
be calculated from this reference pressure. Then 93 can be calculated from Equation (3), 
and T;,,, and T3,,; can be determined from Equations (2) and (4). The other two-thirds of 


the representative cycle are solved in a similar manner. 
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Ist order ode 


Equation (6) 
2nd order ode 










Equation (2) (120-240) 
Ist order ode 








Equation (4) (240-360) 
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Figure 12. Torque Prediction Flowchart 
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Figure 13. Cylinder Gas Pressures (1000 RPM, 100 Ft*Ibf) 


Although the ports are not actually uncovered for the first and last 10 degrees, 
measurements show the assumption of reference pressure is reasonable for the entire 120 
degrees. The assumptions used here limit the feasibility of this solution method to two- 
stroke engines with three or fewer cylinders. 

b. Four-Stroke Engines 

For a four stroke engine, two full rotations of the crankshaft must be 
considered to cover the power, exhaust, intake, and compression strokes. Over the 
representative two rotation (one cycle) period, a cylinder’s pressure can be assumed to be 
equal to intake manifold pressure during the intake stroke, and exhaust back pressure 


during exhaust stroke. Therefore, for a particular cylinder, torque is known for half of the 


Pag 


cycle. The solution method above is feasible for four-stroke engines with four or fewer 
cylinders. A four cylinder engine would require a seven degree of freedom model, which 
could be explicitly solved as above, using the assumptions for intake and exhaust stroke 
pressure. 
c Additional Assumptions for Multiple Cylinder Engines 
Engines with more than four cylinders can be analyzed by this method if 
further assumptions are made. For instance, the cylinder compression stroke can be 
estimated as a polytropic compression of an ideal gas. For a large engine with cylinders 
that fire simultaneously, the two cylinders could be lumped and considered as one inertial 
mass. Also, a measurement device may be placed internal to the engine to measure 
angular velocity at a third degree-of-freedom. 
Zi Interpolation of Data 
The rotational speed data that are collected at the optical encoder and the flywheel 
consists of information that 1s uniformly spaced in the angular position domain, not in the 
time domain. This arises because of the nature of the data collection (See Section I.C 
and JII.D). The raw data for 0,(t) and @5(t) are converted to values which are evenly 
spaced in the time domain for further numerical analysis (specifically, this 1s useful for 
filtering; see section IV.B.3). This requires interpolation of the raw data. Interpolation is 
accomplished numerically using a cubic spline method. This means that the curve is 
assumed to be a 3™ order polynomial with continuous slope at each of the data points. 


Interpolated values of 0(t) are determined for a selected evenly-spaced time basis. 
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5: Signal Filtering by Fast Fourier Transform 





The measured data for position, 6(t), and angular velocity, , contain high 


d6(t) 
t 


frequency components (where “high frequency” refers to frequencies much higher than 
about three times the rotational speed). These frequency components are due to high 
frequency torsional vibration of the crankshaft at the various natural frequencies, and 
random noise from unknown sources. For solution of the equations necessary to predict 
torques, the position data must be differentiated once to determine angular velocity and | 
twice to determine angular acceleration (see section [V.B.1). If raw data were used in the 
analysis, the high frequency components would be greatly amplified by subsequent 
differentiation. However, the torques of primary interest in this problem oscillate at 
about three times the rotational velocity of the engine. Therefore, the much higher 
frequency components are filtered out before differentiation in order to increase the 
signal-to-noise ratio in solving the equations for torque. 

A phase deviation signal is derived from the raw data by comparing it to the mean 
rotational speed. This phase deviation signal is then filtered numerically using a Fast 
Fourier Transform (FFT), removing the high frequency data, and then performing an 
inverse FFT to achieve the desired filtered phase deviation signal. The cut-off frequency 
typically used was between six and nine times the rotational speed of the engine. 

4. Numerical Differentiation 

Once the angular position data has been interpolated and filtered, it must be 
differentiated once to obtain angular velocity and twice to obtain angular acceleration as 
functions of time. A central difference technique is used, where the numerical derivative 


at a point is the sum of the two adjacent differences divided by twice the time difference. 
Zo 


The endpoints are exceptions; they are determined by the single adjacent difference 
divided by the time difference (essentially a forward difference for the first point and a 
backward difference for the last point). No further filtering is required after 


differentiation. 


C. TEST OF SOLUTION METHODS 


The methods previously discussed are used to solve the equations of motion in 
both directions. Using a set of measured cylinder indicated pressures, the accuracy of the 
numerical methods can be tested. The ttme-marching ODE method was used to solve for 
the crankshaft time-resolved angular positions from the measured cylinder indicated 
torques. The predicted 8; and 05 values were then used in the inverse solution method to 
predict the cylinder indicated torques. Comparison of the resulting predicted torques to 
the original measured torques can be used to quantify the accuracy of the numerical 
solution methods. The test results for individual cylinder gas torques are shown in Figure 
i4, and the results for total gas torque is shown in Figure 15. These results show that the 
numerical methods tend to introduce a 2% peak-to-peak error and a 4 degree lag in the 


predicted total gas torque. 
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Figure 14. Test of Numerical Methods for Individual Torques 
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Figure 15. Test of Numerical Methods for Total Torque 
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V. RESULTS 


AG CALIBRATION OF INERTIAL MODEL 


Although the stiffness and inertial parameters of the equations of motion can be 
analytically estimated, some fine-tuning of the values is required to ensure the model 
matches the actual crankshaft. Given the measured torques in the three cylinders, and the 
methods of the previous chapters, the equations of motion may be solved for 8; through 
Q¢. Experimental data is collected for the two measurement points, 8; and 85. Two 
useful comparisons between the measured and predicted time-resolved angular positions 
are the Phase Deviation Plot and the Crankshaft Twist Plot. 

The Phase Deviation Plot shows the oscillation of the angular position about a 
theoretical mean rotating position. It 1s calculated as follows: 

E(t) = 6(t)-G@t (17) 
where @ is the mean rotational velocity. This is the same comparison plot used by 
Hudson [Ref 18]. For a shaft rotating at a steady angular velocity, the phase deviation 
would be zero. The measured phase deviation shows the crankshaft position advancing 
during the power stroke of each cylinder, then retreating during the subsequent 
compression of the next cylinder. A comparison of the measured and predicted phase 
deviation is shown for the crankshaft nose and the flywheel (Figure 16). Plots for other 
runs are found in Appendix D. The phase deviation plot is particularly sensitive to 
assumed values of friction and load torque, and is also useful for validating the model 
inertias. For the 1000 RPM 100 Ft*lbf data, phase deviation shows a maximum error of 
about 9% at the peaks, which correspond to the cylinder compression strokes. Generally, 


the model predicted phase deviation is within 5% of the measured phase deviation for 
33 


both the flywheel and the crankshaft nose. The errors are calculated as a percentage of 
the maximum variation 1n phase deviation. 

The Crankshaft Twist Plot is simply the value of @)-95 as a function of time. It is 
a measure of the total twist along the entire length of the crankshaft. A comparison of 
measured and predicted crankshaft twist is shown in Figure 17. This plot is particularly 
useful for validating the magnitudes of the torsional rigidities used in the model. The 
peak values of the twist occur at the peaks of the gas torque, during the power stroke for 
each cylinder. As expected, the amount of twist is largest for cylinder #1, the farthest 
from the flywheel, and successively lower for cylinders #2 and #3. Correct values of 
torsional rigidity for the model should result in predicted twists comparable to the 
measured value. Comparison of the model predicted twist and the measured twist is 
shown in Figure 17, with additional plots in Appendix D. A max error of about 17% is 
seen at the peak twist values, corresponding to the cylinder power strokes. Generally, the 
model predicted twist is within about 5% of the measured twist. 

As discussed in Appendix C, analysis of the crankshaft natural frequencies can 
also be used to validate the model parameters. The analytical values derived for the 
crankshaft inertias are considered to be reasonably accurate, so the only parameters to be 
adjusted are the torsional rigidities. These are set by the comparison of measured and 
predicted natural frequencies and modes as discussed in Appendix C. Further arbitrary 
adjustment of the parameters to correct the errors in the Phase Deviation and Crankshaft 


Twist plots is not supported. 
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Figure 16. Phase Deviation (1000 RPM, 100 Ft*Ibf) 
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Figure 17. Crankshaft Twist (1000 RPM, 100 Ft*Ibf) 
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|B COMPARISON OF PREDICTED AND MEASURED TORQUES 


Inverse solution of the equations of motion allows prediction of the cylinder gas 
torques, given the time-resolved angular position at two points, 6; and 95. Comparison of 
measured and predicted gas torques for the individual cylinders 1s shown in Figure 18. 
Comparison is not good for individual torques. Besides quantitative errors of over 50% 
at certain points, the predicted torques show misplaced and inappropriate peaks. 
However, it appears that the errors for a particular predicted cylinder torque have 
corresponding offsetting errors in the predicted torque for the other cylinders. This 1s 
evident when measured and predicted total gas torques are compared (Figure 19). The 
predicted total gas torque shows peak-to-peak errors of less than 5%, plus a phase lag of 
about 5-15 degrees. Some of this error is from the numerical solution methods, as 
discussed in the previous chapter. This agreement is good enough to be used for 
localized fault detection. For this particular run of the engine, cylinder #1 gas pressure Is 
significantly lower than the other two cylinders, and the predicted results detect this 


anomaly. Plots for the other data runs are contained in Appendix D. 
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Figure 18. Individual Cylinder Gas Torques (1000 RPM, 100 Ft*Ibf) 
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Figure 19. Total Gas Torque (1000 RPM, 100 Ft*lbf) 
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VI. SUMMARY, CONCLUSIONS, AND RECOMMENDATIONS 


A. SUMMARY 


A three cylinder, two stroke Diesel engine was instrumented with a proximeter 
and an optical encoder for time-resolved angular position measurement at the flywheel 
and crankshaft nose. A torsional model for the engine crankshaft was developed and the 
corresponding equations of motion were formulated. Two separate numerical solution 
methods were developed to solve for the angular positions, given the measured cylinder 
gas torques. These methods were used to calibrate the parameters of the torsional model. 
An inverse solution method was devised to determine the cylinder gas torques, given the 
time-resolved angular positions at two of the degrees of freedom; the flywheel and the 
crankshaft nose. This inverse solution method was shown to be applicable for two-stroke 
engines of three or fewer cylinders, or for four-stroke engines of four or fewer cylinders. 


The predicted cylinder gas torques were compared to measured cylinder gas torques. 


B. CONCLUSIONS 


The torsional model accurately describes the dynamics of the actual crankshaft. 
Experimental] data demonstrated that the model correctly predicted phase deviation at the 
crankshaft endpoints with an error of less than 5%. The model predicted crankshaft twist 
with an error of less than 20%. Predicted natural frequencies from the model agreed with 
the measured frequency spectrum to within 5 Hz for the three vibration modes observed. 

The Finite Element Method (FEM) for direct integration of the equations of 


motion agreed with the Time-marching ODE method to within 1%, and it reduced 


4] 


computational time by a factor of 100 because no iterations were necessary. It was used 
as the primary means for direct integration of the equations of motion in this study. 

The inverse method for predicting cylinder gas torques showed significant errors 
in predicted individual cylinder gas torques. Quantitative errors of over 50%, as well as 
significant wave shape errors, make this method inadequate for reliable prediction of 
individual cylinder torques. However, it is the author’s opinion that this error originates 
in the numerical method used. Specifically, signal filtering tends to create errors in the 
endpoints for the representative cycle. Further analysis of this problem may result in 
successful prediction of individual cylinder torques. 

The inverse method is successful in predicting total cylinder gas torque. 
Predicted total gas torque errors were less than 5%, with slight phase lag errors of 5-15 
degrees. The predicted total gas torque successfully detected a low pressure in cylinder 


#1, showing that the method 1s capable of localizing certain faults to a particular cylinder. 


c. RECOMMENDATIONS 


The failure of the inverse method to predict individual cylinder torques is most 
likely due to problems with the numerical solution methods.. The FFT signal filtering 
induces some errors, which were not sufficiently corrected. The engine speed was not 
exactly steady during data collection, so the filtering process backs out a monotonic trend 
in the phase deviation. That is, the filtered phase deviation is no longer exactly cyclic. A 
correction of some sort should be made for this linear error. The process of signal 
filtering also tends to alter the endpoints of the representative cycle, inducing significant 
errors in the calculated torques at the endpoints. An alternative signal filtering process, 


which avoids these errors, might correct the errors in the results. 
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The assumption of constant parasitic force to model piston ring friction is not 
validated. The presence of crank-angle specific friction "sticking" points would have a 
significant effect on the results. One solution would be to use a motoring dynamometer 
on the engine to produce a graph of the engine friction as a function of crank angle. Also, 
a more detailed analysis of theoretical piston ring friction could lead to more accurate 
modeling. 

From measurements obtained in this study, there is an unknown fault causing 
cylinder #] to have a lower gas pressure than the other two cylinders. As a first step, the 
fuel injectors for cylinders #1 and #2 should be swapped in order to determine the cause 
of the low pressure in cylinder #1. For follow-on experimentation, data should be 
collected for engine runs with known faults. For instance, a defective fuel injector could 
be installed to test the method’s ability to detect a specific fault. 

The use of angular speed measurement internal to the engine would expand this 
method to engines with more cylinders. Although this would be a difficult process for an 
existing engine, mass-produced engines might have such an internal detector installed for 
relatively little extra cost. 

The method for determining TDC for cylinder #1 is inadequate. Currently, the 
procedure of Appendix F, Ref [18] is used to orient the TDC signal on the encoder to 
TDC for the engine. But TDC for the engine is established by a mark inscribed on the 
crankcase and the forward counterweight on the cam follower shaft. Although this shaft 
is directly geared to the crankshaft, gear backlash results in an error of one or two degrees 


when the engine is rotated to TDC. This small error has a significant effect on the 
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magnitudes of the cylinder indicated torques. A better method would be a mark inscribed 
on the flywheel and shroud to ensure correct TDC alignment. 

A sprocket and proximeter assembly might be a more useful means of collecting 
data at the front end of the crankshaft. A 42 tooth sprocket has been obtained which may 
be mounted on the crankshaft nose with the pulley mounting bolt. Since the number of 
teeth on the flywheel (126) is a whole number multiple of 42, a precise alignment could 
be made to calibrate the static phase difference between the two ends of the crankshaft to 
zero. Then the instantaneous twist of the crankshaft could be very accurately measured 
during operation, similar to the method described by Mauer and Watts [Ref 2]. 
Additionally, this data collection method would eliminate the natural frequency torsional 


vibrations of the optical encoder and flexible coupling. 


44 


APPENDIX A. GEOMETRY OF ROTATING COMPONENTS 


While a mass-elastic model has already been developed for this crankshaft [Ref 
23], it is necessary to independently calculate the model parameters in order to verify 
them, and to account for differences in the specific engine used in the research. This 
analysis was carried out using methods and equations from Wilson [Ref 24]. The 


descriptions and values for certain constants used in subsequent equations are listed in 






















Table 6. 
Table 6. Equation Constants 

|Symbol__| Description | Value 
Dy, | Diameter, journal bearing | 3.00in_ 
g __| Gravitational Acceleration _| 386 in/sec™ | 
WK |sRadius ofsGyrationame. cme wali owe 
DRO ery cee me 

Mee, Length, crankpin a0 hee 
(Rs Crankpineccentricity | 2.25in 
R,__| Counterweight Inner Radius | 1.80in 
/Ro _| Counterweight OuterRadius_ | 4.02in 
Me, (Vii | aaa 
'p ___—'| Specific weight of steel__——_| 0.283 Ibf/in” 

7 Rotating Inertia 


An arbitrary objects mass polar moment of inertia 1s calculated as: 


WK, 
J =— (Al) 
g 


where W is the weight of the object, K,,, is the radius of gyration, and g is acceleration 


due to gravity. 
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The polar moment of inertia for a main journal bearing is determined as for a 


solid circular cylinder whose axis is aligned along the axis of rotation: 


WKkee SAD. 1. 
Jin = aa . oo (A2) 


Crankpins are also determined as a solid circular cylinder, but the cylinder axis is parallel 


to and offset from the axis of rotation by the eccentricity: 


WK? De 
J, =— = me be pe . | (A3) 
g g 


As seen in Figure 20, the Crankwebs are of three distinctive types. Crankwebs of 
type (a) connect crankpins for cylinders 1 and 3 to the outer journal bearings, and 
crankwebs of type (b) connect those same crankpins to the inner journal bearings. 


Crankwebs of type (c) connect the crankpin for cylinder 2 to the inner journal bearings. 





(b) (C) 
Figure 20. Crankweb Forms. From Ref [25] 





For all three types, the core geometric shape is an ellipse with one focus centered on the 
journal bearing and the other focus centered on the crankpin. From Table 10.13 of 


Wilson [Ref 24], the radius of gyration for an ellipse can be calculated as: 
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(a 2 axle ) 
cipse = +c (A4) 


where a and b are the minor and major axis, respectively, and c is the offset of the center 
of gravity from the axis of rotation. J can then be calculated for the elliptical portion 
from equation (Al). The counterweight lobes are then treated as semi-circular segments, 


and their contribution to mass polar moment of inertia 1s: 


{ 


1 =o fe _ (A5) 
§ 


where @ is the angle subtended by the counterweight lobe; 120° for crankweb (a) and 70° 
for crankweb (b). 

Rotating inertia for the dynamometer is taken from Ref [26]. The coupling shaft 
between the flywheel and the dynamometer is calculated as a circular cylinder using 
Equation (Al). 

Zz. Torsional Rigidity 

For modeling of the torsional rigidity, the components of the crankshaft must be 
mathematically converted to an equivalent shaft of a constant diameter. This is a simple 
geometric problem for static twisting of the crankshaft, but becomes very complex when 
considering the dynamic crankshaft twist during engine operation. 

Wilson [Ref 24] presents a derivation of the torsional rigidity for various 


crankshaft components. The rigidity for a solid cylindrical shaft of diameter D 1s: 


ee mD*G 
32L, 





(A6) 


This equation can be used to determine the torsional rigidity of any solid cylindrical 
component, including the journal bearings and crankpins. For the crankweb, an effective 
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length is derived based on the bending theory of beams. The crankweb is treated as a 


beam subjected to a bending force by the torsion on the crankshaft: 


0.942R 
Dirge az D| 3 (A7) 
DW 





The equivalent shafting between concentrated inertia points in the model can be 
determined by adding up the torsional rigidity of their components. However, the 
equations above assume an unconstrained crankshaft deflection due to an applied static 
torque. The true equivalent length of the crankshaft elements will be modified by several 
other factors: local deformation where the journal bearings and crankpins join the 
crankweb, bearing restraint on the journal bearing, and non-ideal lever arm at the 
crankweb because the bearing and crankpin are not attached at a single point. [Ref 24] 

Two empirical relations are considered in this study to account for increased 
rigidity of the crankshaft elements due to dynamic constrained operation while mounted 


in the engine block. The first was devised by Carter [Ref 27]: 


LL, +037 Ve Aetna ie 
Low = pa ees a o' 4 , ae a D'| : (A8) 








jb 


cp 


and the second by Wilson [Ref 24]: 


R—-OO(apaD 
ee =| Se (A9) 
wh'" wh 


ejb 4 


jb cp 


F = oe saree =o oe 
ccp 





Table 7 compares the results obtained by each of the methods so far described with the 
values from the model supplied by Detroit Diesel [Ref 23]. The general assumptions for 
the two empirical methods are clear: Carter’s method assumes a stiffer crankpin and more 


flexible crankwebs while Wilson’s method is reversed [Ref 24]. The determination of 
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torsional rigidities for the model is not straightforward; it will require some calibration 


from measured data. 


Table 7. Calculated Torsional — 


Torsional Rigidity (10° Ibf*in/rad) by method: 
=k eva MMR 38g 
MS ae 





12.15 13.50 


3: Auxiliary Loads 

The auxiliary loads, with the exception of the oil pump, are driven off the timing 
gear, just forward of the flywheel. For the model, the total of the auxiliary load torque is 
considered to be placed at 85. The contribution of the individual loads can be determined 


by deriving the power required, then the torque is related to the power P by the equation: 


P 
T =—— (A10) 
22N 


Two cam shafts are driven off the timing gear, at the rear of the crankshaft and 
just forward of the flywheel. The load due to these camshafts primarily results from three 
components: the bearing friction, the operation of the fuel injector pistons, and the 
compression of springs associated with the injectors and the valves. Bearing friction is 
determined in the next section. The action of the injector pistons is considered as a 
polytropic isothermal process of compression from 50 to 2800 psig. Work required to 
compress the springs 1s calculated from: 


ee 
W.. = (0; —6,) (A11) 


Spring 2 Kee 


49 


where 6; and 062 are the initial and final spring deflection, respectively. The load torque 
for the cam shafts 1s calculated from the work done by the shafts over one rotation. 

The fuel pump 1s a positive displacement gear type pump that provides a flowrate 
Q of 1 gpm at 65 psi when the engine is at 2800 RPM. Assuming a pressure at the input 
of about 15 psi, the torque required can be determined from the power P: 

P=QAp (Al2) 

Since the flowrate Q is proportional to the speed N, the fuel pump torque will be a 
constant value regardless of engine speed. 

The water pump is a centrifugal pump that provides 37 gpm at 2800 RPM. Since 
a pressure drop was not provided from the service manual, its effect is estimated as 
comparable to the other auxiliary components. 

A roots blower provides the scavenging pressure that clears the pistons at the 
bottom of the stroke. This blower 1s rated to provide 338 cfm at 2800 RPM. 

The 01] pump is a rotary style positive displacement pump rated to deliver 15 gpm 


at 2800 RPM. The power can be calculated from 
Apg. 
P=mh, = (op) 28) (A13) 


where hag is the head provided by the pump m is the mass flowrate of the oil, p is the oil 
density, and Ap is the pressure increase. The increase in fluid head due to velocity 


increase 1s neglected. 
4. Friction Losses 
Journal bearing friction 1s accounted for by assuming that the two surfaces are 


completely separated by the lubricating film; that hydrodynamic lubrication 1s dominant. 
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Using the relations provided by Heywood [Ref 28], an equation to calculate the parasitic 


torque due to a bearing’s friction 1s derived: 





Di, ae 7p EN W 
- a Jae epi Dl 
ea ee aa 


=—= (A14) 
W.D, ho ie 


where f is the friction factor, Ws is the bearing load, pf is the oil viscosity, h is the 


average bearing clearance, N is the rotation speed in RPM, and Ly and D, are the bearing 


length and diameter. This can be solved for torque to yield: 


3 
TE 6 
4h 


(A15) 
The friction torque is proportional to rotation speed, and independent of bearing load 
under the assumption of hydrodynamic lubrication. This equation can be applied to the 
crankshaft main bearings, the crankpins, and the camshaft bearings. There are four main 
bearings, the coefficients C,, C2, and C3 in the equations of motion are determined as 
one-third of the total friction for the main bearings. Crankpin friction is lumped in with 
the piston ring friction, and the friction due to the camshaft bearings is included with the 
auxiliary loads. 

Piston ring friction is not easily modeled analytically, but is instead estimated by 
empirical methods. From Heywood [Ref 28], piston friction 1s considered as the sum of 
two components: a boundary friction generally proportional to engine loading, and a 
hydrodynamic friction proportional to piston speed. The exact relation is not only 
difficult to predict for any one engine, it will change as the engine’s condition varies. A 
more detailed study of the ring friction 1s beyond the scope of this study. Although it is 


expected that the amount of piston friction varies as a function of the piston speed, a 


constant parasitic force 1s assumed, which corresponds to a parasitic torque Tpar which 


5] 


varies with crank angle, as in Equation (Al6). The derivation of the force/torque 
relationship is detailed in the next Appendix. Since the magnitude of Tpa is small 


relative to T.y, this is considered a reasonable approximation. 


R sin(@ + d) 
cos@ 


(A16) 








The estimates for loads and frictions formulated above are meant to provide a 
relative relation between them. For this study, actual losses were determined from the 


measured pressure data (Table 8). The values for Tioag and Tpar were each estimated as a 


Total Efficiency 
Losses 


C0 
000 [100 [32.4 [587 [2.5__[ 1536 [53.6 65.1% 
150013515371 790 [8242151 ~~ 801 ~~ 
[2000 [160-727 [1033 [1026 [2786 [186 [374% 


fraction of the total losses. 


Table 8. Empirical Friction Losses 


Engine | Dyno Cyl #2] Cyl 
Speed | Load Average | Average 
(ft*lbf) 
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APPENDIX B. GEOMETRY OF RECIPROCATING COMPONENTS 


The nonlinear motion of the piston and connecting rod (Figure 21) presents 
unique geometrical and mathematical problems when modeling a reciprocating engine. 
The following derivations are common in the literature, but are presented here to 


document the detailed analysis required to formulate the torsional model. 





Figure 21. Piston and Connecting Rod. From Ref [19] 


1. Indicated Cylinder Torque 

There are several ways to derive the relation between the gas pressure present in 
the cylinder and the resulting indicated torque applied to the crankshaft. Piston ring 
friction is ignored for this derivation; it is corrected for by a parasitic force as described 
in the previous Appendix. Figure 22 shows the relation of the piston to the crankshaft. A 


Static analysis assumes two forces applied at the piston pin. F, is the net force due to the 
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indicated cylinder pressure, minus the reference pressure applied to the underside of the 
piston from the crankcase (Fp = PrerAp). Additionally, a reaction force F, is applied by the 
cylinder walls on the piston rings; this force constrains the piston to linear motion in the 
cylinder. The resultant force in the direction of the connecting rod is: 


Cc 


i 


cr 


Bl 
CcOsS@ Sate 


where Ap is the cross-sectional area of the piston, P.y is the indicated cylinder pressure, 
and Prep is the reference pressure. The crank angle @ and the connecting rod angles o and 
Y are related by 

Lsin(o) = Rsin(@) and sin(y) = sin(0+0) (B2) 
The resultant torque applied at the crankshaft is then calculated by taking the cross- 
product of the connecting rod resultant force vector and the crankpin position vector, and 
simplifying: 


= 


z NR sin(6 +d) 


COS @ 


(B3) 











R sin y =A, (P., he 





Figure 22. Geometry of Reciprocating Components 
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The same result is also achieved by equating the work done at the piston with the 


work done by the torque on the crankshaft, as demonstrated by Taylor [Ref 29]: 


ii = T’.,40 = a = Pry t (B4) 
Taking the relation for s from Wilson [Ref 24], 
Z oe 
s = Rcos6 + Leos¢ = Rcos6 +|L — Ko sin Ae (BS) 


then 


2.8 ye 
s= 6 —Rsing 7596088 |_ gf _peing—R'sindcos@] og, 
L-—R‘sin’@ cos@ 


and using trigonometric relations, equation B4 will simplify to equation B3. 

2: Reciprocating Torque 

While the rotating parts of the crankshaft maintain a nearly constant angular 
velocity, the reciprocating components are alternately accelerated and decelerated in a 
constrained linear motion. At the crankshaft, this will be seen as a load torque while the 
piston is accelerated from TDC to its maximum speed, and will supply a torque as the 
piston is decelerated to BDC. Taylor [Ref 29] formulates a method of deriving this 
reciprocating torque. | 

First, it 1s necessary to determine the amount of the reciprocating mass. Clearly, 
the entire mass of the piston contributes, but only a portion of the connecting rod 1s 
reciprocating, and the rest must be considered rotating mass. A first approximation 
idealizes the connecting rod as two lumped masses connected by a massless shaft, with 
the same total mass and center of gravity as the real connecting rod (Figure 23). The 
center of gravity for the real connecting rod 1s simply found by balancing the rod; for this 


engine h = 3.5 in. and j = 5.3 in. The portion labeled W, is then added into the 
aD 


reciprocating mass, while the remainder is considered part of the crankshaft rotating 
mass. Then, the instantaneous torque Tye. is found by equating the change in the 


reciprocating mass kinetic energy with the work done by the crankshaft: 


(Fe : = eae 7. i — ™- le, “a 
§ 








W, 


We 


Figure 23. Idealized Connecting Rod. From Ref [29] 


A small correction must then be made to account for the difference between the 
polar moment of inertia for the idealized connecting rod and the actual polar moment of 


inertia: 


W.. |. Rcos@ 

Tree cnr =| Sep “hi P—— (B8) 
| g Leos@ 

where J,, is the polar moment of inertia for the actual connecting rod. While 


Taylor derives series relations to state all values as functions of 0, for this study the angle 


56 


and the distance s are calculated directly, then differentiated numerically within the 
program. 

as Reciprocating Inertia 

During engine operation, the reciprocating components contribute to the rotating 
inertia of the crankshaft. However, the magnitude of the reciprocating inertia varies as a 
function of the crank angle 9. The rotating motion of the crankshaft drives a linear 
motion of the piston in the cylinder. When the piston is at TDC, an incremental rotation 
of the crankshaft results in zero linear motion of the piston, while at 90° the same 
incremental rotation of the crankshaft results in maximum linear motion of the piston. 
Therefore, the influence of the piston mass on the inertia seen at the crankshaft will vary 
during crankshaft rotation. Normally, crankshaft inertial models include an average 
value of one-half the maximum reciprocating inertia to account for the reciprocating 
components. However, in this application, we are interested in crank-angle dependent 
values of torque, so we must account for this crank-angle dependent variation of 
reciprocating inertia. The reciprocating inertia can be calculated as a function of the 


reciprocating mass, the eccentricity, and the crank angle: 


WR 
J =—*_(1-cos(20)) (B9) 
2g 


rec 


The value W,., 1s the weight of the reciprocating components, as determined previously. 
Reciprocating inertia is a separate effect from the reciprocating torque already 

described. The changing value of reciprocating inertia accounts for the extra mass that, 

along with the rotating mass, must be accelerated when the crankshaft is accelerated. The 


reciprocating torque accounts for the energy required to accelerate this reciprocating 


ay 


mass from zero to the current rotational speed of the crankshaft, before it is added to the 


rotating mass. 
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APPENDIX C. NATURAL FREQUENCY AND MODAL ANALYSIS 


A modal analysis can be conducted for the crankshaft using the torsional model. 
Fourier analysis of the torsional vibrations measured at the optical encoder will show 
peaks corresponding to the measured natural frequencies. Comparison of the predicted 
natural frequencies from the model to these measured natural frequencies is a powerful 
calibration tool for fine-tuning the model. 

The matrix equation describing the torsional mode of the crankshaft is repeated 
neue: 

[J]0 +(C]6 +[K]9 = [7] ) 
Neglecting the damping effects, the natural frequencies are calculated by: 

[K]-[7 w’ =[o]=> a’ = eig {JT '[K]} | (C1) 
For the given model parameters (Table 1), the results are tabulated in Figure 24. There 
are six modes of natural vibration, with the first mode being the trivial rigid-body 
oscillation, where there is no crankshaft twist. 

A natural vibration component due to the rigidity and inertia of the flexible 
coupling and the optical encoder disk 1s expected. The given values for the unit are: 

Jrotor = 1.45 x 10° kg*m? 


J coupling =3x 10° kg*m° 


and given 
Tse: 
ie (C2) 
Cnn + ven) 


the natural frequency for torsional vibration for the coupling/encoder unit is: 
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w= 1067 Hz 


Not surprisingly, a large response is seen in the measured data at precisely this frequency. 


ie tt] #2 #3 flywheel dyno 


1679 Hz 1128 H2 429.5 H2 344.4 H2 


1911 H2 





Figure 24. Torsional Vibration Modes 


For highest resolution, a Fast Fourier Transform (FFT) is conducted on the 
measured angular velocity at the optical encoder over the entire 11 cycles of collected 
data. Figure 25 shows the frequency spectrum for measured data from the 1000 RPM, 
100 ft*lbf run. Spectrums obtained for other engine speeds and loads are similar. Figure 
26 is an expanded view of the low-frequency portion of the spectrum. A “comb” of 
amplitude spikes are seen, corresponding to harmonics of the engine rotation frequency, 
16.7 Hz. As expected, a large frequency response is seen at about 1060 Hz, 
corresponding to the natural frequency of the encoder/coupling unit. Additional 


responses are seen at about 430 Hz, 1130 Hz, and 1910 Hz, and these agree with three of 
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the expected modes from the torsional model. A peak seen at about 2400 Hz is from an 
unknown source; it 1s not present in frequency spectrums taken at other engine speeds. 
Natural frequencies that are predicted by the model but not seen in the spectrum 
are probably due to low amplitudes at the crankshaft nose. For instance, the expected 
344 Hz frequency is mostly oscillation of the dynamometer rotor with respect to the 
flywheel; the crankshaft oscillation amplitude would be much smaller. As expected, each 
of the predicted modes has a node close to the flywheel because it contains the bulk of 
the system inertia. Measurements taken at the flywheel would be expected to show 


almost no high frequency vibration, and this is seen in the measured data. 





0 Mode 3 . Mode 4 
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'coupling/rotor 
:OSCillation 
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Figure 25. Frequency Spectrum for Measured Angular Velocity at 9; (0-3000 Hz) 
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Figure 26. Low Frequency Spectrum (0-1000 Hz) 
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Flywheel (Theta five) (rad) 


Crankshaft nose (Theta one) (rad) 


APPENDIX D. ADDITIONAL DATA PLOTS 
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Figure 27. Phase Deviation (1000 RPM, 80 Ft*Ibf) 
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Crankshaft nose (Theta one) (rad) 


aan 1000 RPM, 135 ft*Ibf Phase Deviation 
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Figure 28. Phase Deviation (1000 RPM, 135 Ft*Ibf) 
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ier 1500 RPM, 135 ft*lbf Phase Deviation 
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Figure 29. Phase Deviation (1500 RPM, 135 Ft*ibf) 
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1500 RPM, 160 ft*lbf Phase Dewation 





(pes) (aAy eyay 1) jaaymAl 4 





: oat : : : 
oem ef ee 


+ 





0.01 0.015 0.02 0.025 0.03 0.035 0.04 
Time (sec) 


0.005 








(pes) (2u0 eJay]) BSou YeYysyYUuesD 


0.01 0.015 0.02 0.025 0.03 0.035 0.04 
Time (sec) 


0.005 


160 Ft*Ibf) 


3 


Figure 30. Phase Deviation (1500 RPM 
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2000 RPM, 160 ft*lbf Phase Deviation 





(pei) (any eyayt) jaay 





Mal 4 


0.015 0.02 0.025 0.03 
Time (sec) 


0.01 


0.005 


© 





(pes) (auo ejay) asou WeYsyUuesD 


0.015 0.02 0.025 0.03 
Time (sec) 


0.01 


0.005 


, 160 Ft*Ibf) 


Figure 31. Phase Deviation (2000 RPM 


67 


1000 RPM 80 ft*lbf Crankshaft Twist 
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Figure 32. 
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Figure 33. Crankshaft Twist (1000 RPM, 135 Ft*Ibf) 
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Figure 34. Crankshaft Twist (1500 RPM, 135 Ft*Ibf) 
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2000 RPM 160 ft*!bf Crankshaft Twist 
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Figure 36. Crankshaft Twist (2000 RPM, 160 Ft*Ibf) 
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Figure 37. Individual Cylinder Gas Torques (1000 RPM, 80 Ft*lbf) 
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Figure 38. Individual Cylinder Gas Torques (1000 RPM, 135 Ft*bf) 
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Figure 39. Individual Cylinder Gas Torques (1500 RPM, 135 Ft*lbf) 
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Figure 40. Individual Cylinder Gas Torques (1500 RPM, 160 Ft*Ibf) 
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Figure 41. Individual Cylinder Gas Torques (2000 RPM, 160 Ft*Ibf) 
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Figure 42. Total Gas Torque (1000 RPM, 80 Ft*Ibf) 
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Figure 43. Total Gas Torque (1000 RPM 
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Figure 44. Total Gas Torque (1500 RPM, 135 Ft*lbf) 
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1500 RPM, 160 ft*lbf Total Cylinder Gas Torque 


a =| 








> a Sa aie oe v ae Bs. ar ae i a = 
r [ " ie 
: : ‘ a 8 
‘ : A 
3 — 
é : me ; —* - 
: : om 5 
: e — 
: o pont 
3 neo ee 
: aoe ae 
; of " ¢ 
: ae” ‘ : 
7 caeeet i 
at ee yen * 
ane awe @ee et 
— : . 
oe ow? 
ae eo oe 
ae i : 
- ° a * @ Ravens e+ ; 
eee sear Hn ee eh Sees 
nt ine 8 
‘ero tame : : Sere 
ae ——} me “ 
; : 7 
awe i, 
: oo 
J ae 
‘ =O 
G eae a 
: ese 
: eae el 
, e i 
: “ ww See : 
i » aoe 
: wor 
Se OY 
ae eenee? gene” 
— o 
= weed? 
ov? : 
a ee ee eS as oe; 
—~———e + e 
: = mr a 
'> | 
ee 
oe ees 
eae __f--— 
a sa. 
5 a it = 
neo” 
= ® 
ae ; 
ee? ; 
epee a) = lie 
© © 
ql 
' 


(Jq}.4) anbioy Ses JapUIjAD je}O 1 


/ 
100 150 200 290 300 350 
Crank Angle (deg) 


90 


-400 — 
-600 
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APPENDIX E. MATLAB CODES 


Program used to calculate cylinder gas torques from measured pressures. 


Jo TORQUE36 

Jo This code computes the individual torque contribution of each 
% individual cylinder referenced to TDC of Nr.1 Cylinder. 

J Gas torque is calculated based on measured pressures 


load walk100.md % ECA cylinder #1 pressure data 
load wb1k100.md % ECA cylinder #2 pressure data 
load wclk100.md % ECA cylinder #3 pressure data 


pa = reshape (walk100,5,720); Picyl = [pa(1,:) pa(1)); 
pb = reshape (wb1k100,5,720); P2cyl = [pb(1,:) pb(1)]; 
pe = reshape (wc1k100,5,720); P3cyl = [pe(1,:) pe(1)]; 


eG: % Reciprocating weight (Ibf) 
R22: % Crankshaft Eccentricity (in) 
Bee oor o. % Cylinder Bore (in) 


L = 8.8; % Connecting Rod length (in) 
g = 386; % Gravitational acceleration (Ibf*in/sec*2) 


theta = linspace(0,2*pi,721); % Crank angle vector 
N= 1022; Jo RPM 

Load = 100; % Ft*lbf 

omega = 2*p1*N/60; % rad/sec 


dt = (60/N)/720; 

sl] = R*cos(theta) + sqrt(L’%2 - (R*2)*sin(theta).*2); 

s2 = R*cos(theta-4*pi/3) + sqrt(L*2 - (R*2)*sin(theta-4*pi/3).%2); 

s3 = R*cos(theta-2*pi/3) + sqrt(L*2 - (R“2)*sin(theta-2*pi/3).%2); 

Spl = deriv(s},dt); Spdl = deriv(Sp1,dt);% Piston speed (in/sec) and 

Sp2 = deriv(s2,dt); Spd2 = deriv(Sp2,dt);% Piston acceleration (in/sec’2) 

Sp3 = deriv(s3,dt); Spd3 = deriv(Sp3,dt); 

Tlrec = -(W/g)*Spd1.*(Sp1/omega)/12; 

T2rec = -(W/g)*Spd2.*(Sp2/omega)/1 2; 

T3rec = -(W/g)*Spd3.*(Sp3/omega)/12; 

Fpar = 201; Jo \bf Parasitic force 

Tparl = Fpar*abs(Sp1/omega)/12; % ft*lbf Parasitic torque 

Tpar2 = Fpar*abs(Sp2/omega)/12; 

Tpar3 = Fpar*abs(Sp3/omega)/12; 

pref = 14.706 + 0.00205*(N-634); % psia 

Ticyl=-((Plcyl.*pref-pref) *(p1*B%2/4).*Sp 1/omega)/12; Tool eior 
T2cyl=-((P2cyl.*pref-pref)*(pi*B%2/4).*Sp2/omega)/12; Hose lA 6) 
T3cyl=-((P3cyl.*pref-pref) *(pi*B%2/4).*Sp3/omega)/12; fo) ELBE 

Tbrg = (0.0003403)*N; % FT-LBF Bearing friction per cylinder 
‘Faux = 38; % FT-LBF Valve train and auxilltaries 
Tpmp = 0.792; % FT-LBF Oil pump load 

Ttot = Ticyl+T2cyl+T3cyl-Tpar!-Tpar2-Tpar3-3*Tbrg-Taux-Tpmp+T Irec+T2rec+T3rec;% FT-LBF 
figure (1) 

degrees = theta*1 80/pi; 

plotidecrees Pilcy!. kx deorces,P2cyi, ko degrees, Pacyl. k. ) 

axis([0,360,0,90]) 

title (Cylinder Pressures referenced to TDC of #1 Cylinder’) 

xlabel (Crank Angle (degrees)’) 

ylabel (‘Cylinder Indicated Pressure (bars, absolute)’ 

legend (Cyl#1’,Cyl#2’Cyl#3) 
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grid 

orient landscape 

figure(2) 

plot(theta,T1lcyl, bX’,theta,T2cyl, bO’*theta,T3cyl, b.) 

grid, hold on 

plot(theta,T lrec, rX’,theta, T2rec, rO’*,theta, T3rec,r.) 

plot (theta, Ttot,’g’,theta, Load *ones(size(theta)), ‘c’) 

title(Individual Cylinder Torque input Referenced to TDC of Nr. | Cylinder’? 
ylabel(Gas Torque (FT-LB)’ 

xlabel( Degrees after TDC of NR. | Cylinder’? 

legend(Cyl #1’,’Cyl #2’, Cyl #3) 

orient landscape 

Jo COMPUTE AVERAGE TORQUE CONTRIBUTION OF EACH CYLINDER 
Tlcyl_avg=trapz(theta,T 1 cyl)/(2*pi) 

T2cyl_avg=trapz(theta,T 2cyl)/(2*p1) 

T3cyl_avg=trapz(theta, T3cyl)/(2*p1) 

Avg_Torque_input=trapz(theta, Ttot)/(2*p1) 

Torque_out_to_Torque_in = Load/Avg_Torque_input 


Programs used to calculate angular positions, given cylinder gas torques 


(time-marching method). 


Jo MEASPRD3 

Jo Concurrently plots measured and predicted responses 

%o using time-marching direct integration method 

Jo 

%o Section One plots the measured response 

Jo 

load w11k100.md % Load optical encoder data file 

t= wllk100(C,1); 

tt = diff(t); % determine dt’s 

it = (resnanettt./20/11)). % Phase lock one cycle 

{ttt = mean(tt, 1); % Ensemble average the phases 

position = linspace(0,(2*p1),721); % The known positions of the O.E. windows 
pos = position(1:720); 

omega = (1/sum(tttt))*2*pi; % Mean rotational velocity (rad/sec) 
timem = [0 cumsum(tttt)]; % Time vector corresponding to position 
angposm = position-timem.*omega; % Angular position (radians) 
plot(timem,angposm, rO’) % Plot measured angular position vs time 
hold on 

% 

%o Section Two plots the predicted response based on pressure data 

%o Uses a SIX SECOND ORDER SUMULTANEOUS EQUATION ODE SOLVER 
Jo Calls "deqns.m" which defines the system of 2nd order ode’s 

%o 

global Tload Tlcyl T2cyl T3cyl Taux Tpmp j2rec j3rec j4rec 1 

load walk100.md % ECA cylinder #1 pressure data 

load wb1k100.md % ECA cylinder #2 pressure data 

load wclk100.md % ECA cylinder #3 pressure data 


pa = reshape (walk100,5,720); pa = pa(1,:); 
pb = reshape (wb1k100,5,720); pb = pb(1,:); 
pe = reshape (wc 1k100,5,720); pc = pc(1,:); 
Plcyl = [pa(1,:) pa(1)];% Cylinder pressures'(bars) 
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P2cyl = [pb(1,:) pb(1)]; 

P3cyl = [pce(1,:) pe(t)]; 

shp = [00000 -1.8}*le-03; 
Sveti lel 1106.3: 

ic = [shp shv]; 

Tload=1200; 
timep=linspace(0,sum(tttt), 721); 
angposp = zeros(1,/721); 


% define IC’s 

% \bf*in Load torque 

% divide one rev into 720 divisions 
% initialize predicted position vector 


Cylinder bore 

Crankshaft eccentricity 

Weight of reciprocating components 
gravitational acceleration 
Connecting rod length 


B=3.673. % in 

R= 225: % in 

W = 72556: % \bf 

c= 356: % in/sec*2 
[= 75758 % in 

N = omega*60/(2*p1); % RPM 


pos = linspace (0,(2*p1),721); 


pref = 14.706 + 0.00205*(N-634); % psia 


dt = (2*pi/omega)/720; 


sl = R*cos(pos) + sqrt(L*2 - (R%2)*sin(pos).42); 
s2 = R*cos(pos-4*pi/3) + sqrt(L*2 - (R*2)*sin(pos-4*pi/3).%2); 
s3 = R*cos(pos-2*pi/3) + sqrt(L42 - (R42)*sin(pos-2*pi/3).%2); 


Sp! =derivigidt): Spdl = derivGaldn: 
Sp2 = deriv(s2,dt); Spd2 = deriv(Sp2,dt); 
Sp — deriviss di): sopds = deningsp3.dt); 
Fpar = 100; % 
Tparl = Fpar*abs(Sp!/omega); %o 
Tpar2 = Fpar*abs(Sp2/omega); 

Tpar3 = Fpar*abs(Sp3/ome ga); 

pos = linspace (0,(2*p1),721); 


pref = 14.706 + 0.00205*(N-634); % psia 
%o\bf*in Oil pump torque 


Tpmp = 9.5; 


% Piston speed (in/sec) and 
% Piston acceleration (in/sec*2) 


lbf*in 
Ibf*in 


Parasitic force 
Parasitic torque 


Tlcyl = -(Plcyl.*pref-pref)*(pi*B%2/4).*Spl/omega; % (in*lbf) 


T2cyl 


-(P2cyl.*pre f-pref)*(pi*B%2/4).*Sp2/omega; % (in*Ibf) 


T3cyl = -(P3cyl.*pref-pref)*(pi*B%2/4).*Sp3/omega; % (in*Ibf) 


Taux = 1608: 
j2rec = (W*R42/(2*g))*(1-cos(2*pos)); 


j3rec = (W*R‘2/(2*g))*(1-cos(2*(pos-4 *pi/3))); 
j4rec = (W*R42/(2*g))*(1-cos(2*(pos-2* pi/3))); 


Tlrec = -(W/g)*Spd1.*(Spl/omega); 
T2rec = -(W/g)*Spd2.*(Sp2/ome ga); 
T3rec = -(W/g)*Spd3.*(Sp3/omega); 
Tlcyl = Tlcyl - Tparl + Tlrec; 
T2cyl= P2eyl- Tpar2 + P2rec; 
Tey W3ey= Tpar3"+ 1 3rec; 

Step = 1; 

theta = zeros(721 ,6); 

theta(1,:) = shp; 

[OL | step: 720: 


%lbf*in Valvetrain and auxilliary torque 


Jolbf*in*sec*2 
%olbf*in*sec*2 
%olbf*in*sec42 


(T,x] = ode45(deqns’[timep(i) timep(i+step)],ic); 


ic = x(length(T),:); 
theta(i+step,:) = x(length(T), 1:6); 


% reinitialize IC’s from previous iteration 


angposp(i+step) = x(length(T),1) - omega*timep(itstep); 


end 


plot(timep,angposp, bX’) 9% Plot predicted angular position vs time 
title( Time Marching ODE45 Method Comparison to Measured Data’) 


xlabel(’Time (sec)? 
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ylabel( Theta one Phase Deviation (radians)? 
legend (Meas ’,Pred =’ 

grid 

orient landscape 

shaft_pos = (ic(1:6)-2* pi)’ 

cycle_stat = (ic - [(shp+2*pi) shv])’ 

% Plot predicted angular velocity 

figure (2) 

omegap = [shv(1) (diff(angposp)./timep(2)+omega)]; 
omegam = (1./(/20*tttt))*2*pi; 

omegam = [omegam omegam(720)]; 
plot(timep,omegap, b’,timep,omegam, T’7) 
hold on 
plot(timep,omega*ones(size(timep))) 
xlabel(’time(sec)’) 

ylabel(’Angular Velocity (rad/sec)) 

grid, orient tall 

figure (3) 

plot(timem,angposm, k.) 

hold on 

plot(timep,angposp, k’) 

title’ Time Marching ODE45 Method Comparison to Measured Data’) 
xlabel( Time (sec)? 

ylabel( Theta one Phase Deviation (radians)? 
legend (Meas °Pred 

grid 


% DEQNS function to determine six second order differential equations 
% To be used in ode45 fctn in measpred program 


Ty se a ne 
function xdot=deqns(t,x) 

Oe a cncictinn 

global Tload Tlcyl T2cyl T3cyl Taux Tpmp j2rec j3rec j4rec 1 
Jo------ Constants to be used in differential equations------ 
j1=.02443; %\b*in*sec*2/rad 

j2=.2482; %|b*in*sec*?2/rad 

jo= 1462, Jolb*in*sec*2/rad 

j4=.2482; Jolb*in*sec*2/rad 

jo=).222. Jolb*in*sec’2/rad 

}6=Zsr0: %\|b*in*sec*2/rad 

kl=3 lke6; Jolb *in/rad 

k2=7.00e6; %\b*in/rad 

k3=7.00e6; %\b*in/rad 

k4=10.82e6; Jo\lb*in/rad 

k5=1.304e6; J |b*in/rad 

C= 01 J%|\b*in*sec/rad 

C23=01 %\b*in*sec/rad 

e34=01- % |b*in*sec/rad 

e45— 017 %\b*in*sec/rad 

e56=.01; J |b*in*sec/rad 

e2=0:013. %|b*in*sec/rad 

e3=0 Ole; %|b*in*sec/rad 

c4=0.073; %\b*in*sec/rad 

Of yaaa es ca eee eg 

% 12 first order equations which define the 
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Jo original 6 second order equations of motion 


xdoUl)=x(/); 

XGdo0(2)=x(3); 

xdo(3)=x(9); 

xdo(4)=x(10); 

xdo(5)=x(11); 

xdo(6)=x(12); 

xdo(7)=((c12*x(8))-(c 12*x(7))+(k1 *x(2))-(k1 *x(1))-Tpmp)/j1; 

xdo(8)=(T I cyl(i)+(c23 *x(9))-((c23 +c 12+c2)*x(8))+(c 12*x(7))+(k2*x(3))... 
-((k1+k2)*x(2))+(k1 *x(1))V/Q24j2rec(1)); 

xdo(9)=(T2cyl(i)+(c34*x(10))-((c344+c234+c3)*x(9))+(c23*x(8))+(K3 *x(4))... 
-((k2+k3)*x(3))+(k2*x(2)))/Q34j3rec(1)); 

xdo(10)=(T3cyl(i)+(c45*x(11))-((c45+c34+c4)*x(10))+(c34*x(9))+(k4*x(5))... 
-((k3+k4)*x(4))+(k3*x(3)))/Q4+4j4rec(1)); 

xdo(11)=((c56*x(12))+((c56+c45)*x(11))+(c45*x(10))+(k5 *x(6))... 
-((k4+k5)*x(5))+(k4*x(4))-Taux)/j5; 

xdo(12)=(-Tload-(c56*x(12))+(c56*x(11))-(k5 *x(6))+(k5 *x(5)))/j6; 

xdot=xdo’, % vector defining the equations of motion 


Programs used to calculate angular positions, given measured cylinder gas 


torques (finite element method). 


Jo MEASES 

% Determines and plots comparisons of the following: 

% (1) Measured response from flywheel and optical encoder data 

% (2) Predicted response from inertial model based on measured 

% pressure data from the three cylinders 

% Evaluates predicted response using a finite element formulation 

Yon--2-------------------------------------- 

Jo Section One plots the measured response 

CMe a eRe ok hoe 

%-----Measured response from flywheel----- 

load w51k100.md % Load time data from MDA 

t = [0;w51k100(1:7938,1)]; % Extract me data only 
_tt=difi(t); % COMPUTES THE DEL_T’S 

tt = (reshape(tt, 126,63))’; % Phase lock one cycle 

tttt = mean(tt); % Ensemble average the phases 


dtrat = tt(1)/mean(tt(2:63,1)); 

teeth = linspace(0,2*pi, 127); tooth = 2*pi/127; 

th5 1 = tooth*dtrat; 

pos) = [0 teeth(1:125)+th5 1]; 

pos5(127) = 2*pi; 

time5 = [0,cumsum<(tttt)-(1-dtrat)*tttt(1)]; 

time5(127) = sum(tttt); 

omega = (1/max(time5))*2*pi1; 

angpos) = pos5-time5.*omega;%Angular position (radians) 


figure (1) 

subplot (2,1,1) 

plot(time5,angpos5, k.’) % Plot measured angular position vs time 
hold on 

%-----Measured response from optical encoder----- 

load w11k100.md % Load time data from MDA 

C= (0-1 KO0Cl 7920°1)]. % Extract time data only 
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tt = diff(t); % COMPURES THE DEES TS 


it = (reshape(tt, 720M) % Phase lock one cycle 

{ttt = mean(tt); % Ensemble average the phases 

pos1 = linspace(0,(2*p1),721); % Crank angle position for O.E. windows 
time1l = [0 cumsum(tttt)]; % Time vector corresponding to position 

time] = timel *(max(time5)/max(timel)); | % Adjust times to agree 

angposl = pos1-timel.*omega; % Angular position (radians) 

figure (1) 

subplot (2,1,2) 

plot(time1,angpos1,k.’) % Plot measured angular position vs time 
hold on 

Ole ie ook oo fom eee si bess re... ....... 

Jo Section Two plots the predicted response based on pressure data 

Oj, oe ee Me eid 25 os. 


global j2rec j3rec j4rec 
Jo------ Cylinder pressures (bars) ------ 


load walk100.md % ECA cylinder #1 pressure data 
load wb1k100.md % ECA cylinder #2 pressure data 
load wclk100.md % ECA cylinder #3 pressure data 


pa = reshape (walk100,5,720); Plcyl = [pa(1,:) pa(1)}; 

pb = reshape (wb1k100,5,720); P2cyl = [pb(1,:) pb(1)]; 

pc = reshape (wce1k100,5,720); P3cyl = [pc(1,:) pc(1)]; 

% Nariable descriptions 

Yo k=element matrix 

% f=element vector 

% kk =compressed system matrix 

% ff =system vector 

% bcdof =a vector containing dofs associated with boundary conditions 
% bcval =a vector containing boundary condition values associated with 
Jo the dofs in bcdof’ 


ee 
% input data for control parameters 

Jo ----- === nn nnn nnn nnn nnn === 

neli= 720: % number of elements 

nnel=_2. % number of nodes per element 

ndof = 6; % number of dofs per node 

nnodem 721; % total number of nodes in system 

sdof = nnode*ndof; % total system dofs 

Tload = 1200; % (Ibf*in) 

R=) J (in) Crankshaft eccentricity 

W =7.556; % (\bf) Weight of reciprocating components 
g = 386; % (in/sec’*2) Gravitational acceleration 

B= 3.675: % (in) Cylinder Bore 

E=s 3. % (in) Connecting Rod length 

N = omega*60/(2*pi); J RPM 


dt = (2*pi/omega)/720; 

sl = R*cos(pos1) + sqrt(L%2 - (R%2)*sin(pos1).%2); 

s2 = R*cos(pos1-4*pi/3) + sqrt(L’2 - (R%2)*sin(pos 1-4*pi/3).%2); 

s3 = R*cos(pos1-2*pi/3) + sqrt(L’*2 - (R*2)*sin(pos1-2*pi/3).%2); 

Spl = deriv(sl,dt); Spdl =deriv(Spl,dt);  % Piston speed (in/sec) and 
Sp2 = deriv(s2,dt); Spd2 =deriv(Sp2,dt); | % Piston acceleration (in/sec’2) 
Sp3 = deriv(s3,dt); Spd3 = deriv(Sp3.dt); 

Fpar = 100; Jo Ibf*in Parasitic force 

Tparl = Fpar*abs(Sp1/omega); Jo Ibf*in Parasitic torque 

Tpar2 = Fpar*abs(Sp2/omega); 
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Tpar3 = Fpar*abs(Sp3/omega); 

pos = linspace (0,(2*p1),721); 

pref = 14.706 + 0.00205*(N-634); % psia 
Tpmp = 9.5; %olbf*in Oil pump torque 


Tlcyl = -(Plcyl.*pref-pref)*(pi*B%2/4).*Sp1l/omega; % (in*lbf) 
T2cyl = -(P2cyl.*pref-pref)*(p1*B%2/4).*Sp2/omega; % (in*lbf) 
T3cyl = -(P3cyl.*pref-pref)*(pi*B%2/4).*Sp3/omega; % (in*lbf) 
Taux = 168; %olbf*in Valvetrain and auxilliary torque 

qetee = (Wk (2-2) { leeos(2 pos): Tolbr*in* sec 2 
j3rec = (W*R42/(2*g))*(1-cos(2*(pos-4*p1/3))); %olbf*in*sec*2 
j4rec = (W*R42/(2*g)) *(1-cos(2*(pos-2*pi/3))); Zo\lbf*in *sec2 


Tlrec = -(W/g)*Spd1.*(Sp 1/omega); 
T2rec = -(W/g)*Spd2.*(Sp2/omega); 
T3rec = -(W/g)*Spd3.*(Sp3/omega); 
Tlcyl = Tlcyl - Tparl + Tlrec; 
T2eyl'= T2cyl=fpar2 + T2rec; 
T3cyl = T3cyl - Tpar3 + T3rec; 


pack 

0 

% input data for nodal coordinate values 
CE. no eee ec cekaten das ee 
tcoord = linspace(O,max(timeS),721); 

COE oat AOI ENE 2 LE 
% input data for nodal connectivity for each element 
oe a as ee Os ee Bes se 
nodes = [(1:nel)’,(2:nnode)]; 

Cf, eee iene ee Tee ee 

% input data for boundary conditions 
ON a a et 

Jo Dirichlet Boundary Conditions 


shp = [0 0 000 -200)*1le-05; 
bedof = [1,2,3,4,5,6,sdof-5,sdof-4,sdof-3,sdof-2,sdof- 1 ,sdof]; 
beval = [shp (shp+2*p1)]; 


A hy ee ee ee a a 
% initialization of matrices and vectors 

id eR rae 

ff = zeros(sdof, 1); % initialization of system force vector 

kk = zeros(sdof,18); % initialization of compressed system matrix 
index = zeros(sdof,1); % initialization of kk index vector 

9 Resta Ne | rr a ee ee 


for i=1:nel; % loop for the total number of elements 
nl=nodes(i,1); nr=nodes(i,2); % extract nodes for (iel)-th element 
tl=tcoord(nl); tr=tcoord(nr); % extract nodal coord values for the element 
k = Kelm(tl,tr,1); % compute element matrix 


f = ((tr-tl)/2)*[-Tpmp; Tlicyl(i); T2cyl(i); T3cyl(i); -Taux; -Tload;... 
-Tpmp; Tlicyl(i+1); T2cyl(i+1); T3cyli+1); -Taux; -Tload]; 
% compute element vector 


fOF tal] 2: % assemble element matrices and vectors 
ff((i-1)*6+11) = ff((i-1)*6+i1) + f(ii); 
if 11<=6; 


kk((i-1)*6411,:) = kk((i-1)*6+411,:) + [0,0,0,0,0,0,k(ii,:)); 
index(i*6+11) = (1-1)*6; 
else 
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kk((1-1)*6+n1,:) = kk((i-1)*6+11,:) + [k(ii,:),0,0,0,0,0,0); 


end 

end 
end 
indéx(1:6)'==G67| let 1-1-1]: 
G7, ocianciae cee ee ee 
% apply boundary conditions 
ORs erage ene ds 
% Dirichlet Boundary Conditions 


for 1 = I:length(bcdof); 
kk(bedof(1),:) = zeros(1,18); 
kk(bcdof(i),bcdof(i)-index(bedof(i))) = 1; 
ff(bcdof(1)) = beval(i); 

end 


for 1 = I:(sdof-1); 
f = fix((i-1)/6)*6+ 12; 
if 1>(sdof-6); f=sdof; end 
for 1 = (i+1):f; % other rows with data in column 1 
v = kk(ii,1-index(i1))./kk(i,1-index(1)); % multiplier 
kk(1i,1-index(i1)) = 0; 


for ) = (it+l):f; % data elements in row 1 
kk(ii,J-index(11)) = kk(11,J-index(11))-v*kk(i,j-index(i)); 
end 
ff(ii) = ff(i) - v*ff(i); 
end 
end 
a a eS ee ice 
% Solve matrix eqn for theta 
Jp -nn--nn nanan nena 2-222 -- == == 


theta = zeros(1,sdof); 
theta(sdof) = ff(sdof)/kk(sdof,sdof-index(sdof)); 
for 1 = (sdof-1):(-1):1; 

f = fix((Ci-1)/6)*6+12; 

if 1>(sdof-6); f=sdof; end 
» for) = +1):f; 

ff(i) = ff(i) - kk(i,j-index(1))*theta(j); 

end 

theta(i) = ff()/kk(i,1-index(1)); 
end 
theta = reshape(theta,6,nnode); 
angposp1 = theta(I,:) - omega*tcoord; 
angposp5 = theta(5,:) - omega*tcoord; 


figure (1) 

subplot(2,1,1) 

plot(tcoord,angpospS, k’) 

xlabel(’ Time (sec)’) 

title (1000 RPM, 100 ft*Ibf Phase Deviation’) 
ylabel( Fly wheel (Theta five) (rad)’) 

grid 


90 


subplot(2, 1,2) 

plot(tcoord,angposp1,k’) 

xlabel(’ Time (sec)’) 

ylabel(’Crankshaft nose (Theta one) (rad)’) 


grid 

legend (Meas ’,Pred 

Gt Rees ts 2 re 

% Compare Other degrees of freedom to theta 5 
2 FS. 

twist = theta - ([1;1;1;1;1;1]*theta(5,:)); 

figure (3) 


plot (tcoord,twist(1,:),’y’), hold on 

plot (tcoord,twist(2,:),r’) 

plot (tcoord,twist(3,:),’z) 

plot (tcoord,twist(4,:),b) 

plot (tcoord,twist(6,:),’c’) 

title’ Angular Deviation from Theta Five (Flywheel)’) 
xlabel( Time (sec)’ 

ylabel(’Angular Deviation from Theta Five (radians)’) 
legend (O.E.’,Cyl #1’,'Cyl #2’,’Cyl #3’, Dyno’, Measured’) 
erid 

orient landscape 


thetarl = interp1(time!,pos]1,tcoord,’spline); 
thetar5 = interp1(time5,pos5,tcoord,’spline); 
figure (6) 

plot (tcoord,thetar 1-thetar5, k.’) 

grid, hold on 

plot (tcoord,twist(2,:),k) 

title (1000 RPM 100 ft*lbf Crankshaft Twist’? 
xlabel (Time (sec)’ 

ylabel (Twist (radians)’) 

legend( Meas ’Pred 


function [k] = kelm(tl,tr,1) 

% 

% KELM calculates the element matrix k, 
% given tl and tr as inputs 

% Global variables 

global j2rec j3rec j4rec 


n= tr-tl: 

jl = 0.02443; %\bf*in*sec*2/rad 

SO es %olbf*in*sec*2/rad 

Jo 01462. %lbf*in*sec*2/rad 

j4 = 0.2482; %lbf*in*sec42/rad 

ie ee %|bf*in*sec*2/rad 
jO'= 0.2870; %\bf*in* sec*2/rad 

ki 3 eG; %\bf*in/rad 

kZ= 7 00e6: %\bf*in/rad 

k3 = 7.00e6; %\bf*in/rad 


ot 


k4 = 10.82e6; J%olbf*in/rad 
k5 = 1.30466; %\of*in/rad 


el? = 0.01; %\bf*in*sec/rad 
e23=C01- % |bf* in * sec/rad 
e34=0:015 %\bf*in *sec/rad 
c45 = 0.01; %|\bf*in*sec/rad 
e56:= 0:01 J \bf*in*sec/rad 
eZ =0:013; %\bf*in* sec/rad 
cs = 0012: J \bf*in*sec/rad 
c4 =O 01s %|bf*in*sec/rad 
Oi, eines Sota eee eeaee 

k= Zeros:(Iz,12): 

Of. Saree sea 


KOR D=1/6*(27k I tte*3-2*k) *tl*3-3 *c 1272-37012 *1%2+6* ki irl 2-07 leone 
ql tle67cl2 ste *tl-Ge kl 2h; 

KC 2=NW6 "G2 "Kl * tr 3427 k1 tls 43 tol este*2+3*cl2 *1*2-6*k] *te*tl*2-6*cll2 ie tl 
6* kl ttre 27/2: 

k17 == 1/6*CK1 *tr*3 +k) *1%3-37 12 *tr*2-3*e12 tI 2-3* kl trl ki 2 toa, 
yl ¥tr+6*) 1 *tl+6%c1 2 *tr*tl)/h’%2; 

k(18)==176* (Kk)? tr73-k 1 *tl*343*c 12 *tr*2+3*012 71424 3SRl tel so -3 kl te 2 tl-o7 
eles tre thn: 

KZ, = 1/6" (-2"k 1] str 32 elt 3 + 3 cl 2 tr“ 2 eee 2 * 1)" 2-6*k 1] *tr“tl* 2-67-12 te “tl. 
6* Kite 27th 2: 

k(2,2)=1/6*(6*k2*tr*tl42-3*c2*tr42-3*c2*tl42-6*j2reci) *tr+6*)2rec(i) *tl-6*)2 *tr+... 
6*j2*11+6* tl *tr*¥c23+6*tr*c2 *tl-6 *k2*tl*tr*2+6*%c)] 2 *tr*tl-6*kietr 21142 *.... 
k2*tr®3-2*k2 *1I*3-3*¢23 *tr*2-3 *c237tl*2+2*kl *tr*3-2*K1 see eee licett \2- 
3*C12472+07k) *te* 142 )/h*2: 

KS) 1/67 E222 Fr 32 *k2 114 3-67 k2 ti 3 "e235 ir 2323 tl 2-6F il ih e234 
OPk2 *tl*tr°2)/n*2: 

KOZ. 7 = Ga kil te43-k 1 tl43 43%C12 *te2 243 7c 12 sao 3 tl 2-3 kt olor 
C12 tht ne 2. 

K(2,8)=- 1/6*(-3 *k2 *tr*tl*2-3 *c2 *tr*2-3 *c2 *t142-6* | 2rec( i) tieeom erect atl; j2~... 
tr+6*j2*tl+6*tl*tr*¥c23+6*tr*c2 *tl+3*k2 *tl *tr*24+6¥*c 12 *tr*tl+3*kl *tr%2*... 
Hak? tre 3 4k? stl? ses2e23* ir 2-5 C23 tl) -Kk] tir 34K) 1 3-3 el? tr 2-5... 
Cle tl) 2-3 kl ie tl*2)/h*2: 

KO 9 SG ke tr 3-k 2 114343 *k2 ttr*tl4243462 3 ir 2 +3 7623 711"2-3*k2ttl ir 2 -o 
tht, e232: 

KGa = WW re kK? *tr*3 +2 *k2 * 1143-6 *k2 “tr tl 243 *e23 tr 243 e725 112-6 tl tir c+. 
6*K2 *tl*tre 2/2: 

k(3,3)=1/6*(6Fk2* tr *t142+6 7 tl* tr *623-6*k2*tl*1r*2+6*k3 *tr*1l*2-67k3 *1r*2 11+ oe 
teat 3432 tps ke 1% 3-3 *C23* treo 5 8C2 3h eek tr oak 
11°3=3*634 *tr*2-3*c 347114 2-3%c3 *tr*2-3* eS 11" 2+67()*c3"tr-6*j3"trt0*.. 
j3*tl-6*j3rec(i)*tr+6*j3rec(1)*tl)/h’?2; 

K(G4 = Geek tr 342731183 -67k37tr 719243 4034 81r4243"c34 tl" 2-6 tel eb... 
OF k3s*tr 2h 2: 

KG SS Work te 3-k2 Ft ok? tel 3 "C23 tr 243 e231 2-3 he i 2 
tl*tr*e23) a2: 

kG,9)=- 1/673" k2 * tr * 114246" tl *tr*c2343*k2* fl *te*2-3* k3 tr 11243 ko Ae ato 
tr*th*c34-k2*te*3--k2 *11%3-3 t623 *tr* 2-3 *C2 341192 -kS tir ks tl ec oe 
tr°2-3*¢34*tI*2-3 *c3*tr%2-3*c3 *tl*2+67tl*c3 *tr-6*)3*tr-+6*)3 *tlGame 
j3rec(i)*tr+6*j3rec(1)*tl)/n%2; 

k(3,.10)=21/670K 37 te*3-k3 #04343 *k3 Ftrttl42+3 *e34*tr4 2437634 11h 223 ho 2 lee 
OF tet] 2634 hee. 

k(4,3)=1/6*(-2*k3 *tr43+2*k3*tl43-6*k3*tr*tl42+3 *c34 *tr*24+3*c34 *t142-6* tr l*e34... 
G*k3 *tr 2 thy. 


oe 


k(4,4)=1/6*(6 *tr*c4 *tl+6 *tr*c4 5 *t1+6*k3 *tr*t142-6 *k3 *tr42 *t1+6 *tr*tl*c344+2*k3*... 
tr*3-2*k3*t143-3*c34 *tr42-3 *c34 *t142+6 Ftr*k4 *t1%2-6*k4 * tl *tr42+-2*k4 *... 
tr*3-2*k4 *t143-3*c45 *tr42-3 *c45 *t142-3 *c4 *tr42-3 *c4 * 1142-674 *tr+6*... 
j4*tl-6*)4rec(i)*tr+6*j4rec(i)*tl/h’2; 

k(4,5)=1/6*(-2 *k4 *tr434+2*k4 *t143-6* tr*k4 * 14243 *c45 *tr424+3 *c45 *t142-6* tr*c45*... 
tl+6*k4*tl *tr42 V/h%2; 

KAO) =" 1/6 a 3 es Ss ir I 23 2634774 3 7) 34 1182-3 KS “ir *t1-6*... 
tr*tl*c34)/h‘2; 

k(4,10)=- 1/6*(6* tr*c4 *tl+6*tr*c45 *tl-3*k3 *tr*t]424+3 *k3 *tr*2 *tl+6*tr*tl *c34-k3 *... 
(eS Fo tls 25 e347 2-3 e346 Ul 225th k4 tI 2s ka kK te 
k4*t]43-3*c45 *tr42-3*c45 *tl42-3 *c4 *tr42-3 *c4 *t142-6*)4 *tr+6*j4 *tl-6*... 
j4rec(i)*tr+6*j4rec(i)*t)/h’2; 

KG =- 6k 3-0-3 ka tr 23 24 5 i 437045 0142 +3 tr *k4 #014 2-6*... 
tr*c45*tl)/h*2; 

k(5,4)=1/6*(-2*k4 *tr43 +2 *k4 *t143-6 *tr*k4 *t142+3 *c45 *tr42+3 *c45 *t142-6* tr*c45* tlt... 
6*k4 *tl*tr*2)/h’2; 

KO 53)=1/67 (2 "ka tr 3-2 *k4 Sees tr 3-2 *k5 *t1%3-3* C45 tr 2-3 *c45 * 11% 2-3 *C56*... 
tr92-3* C56 *tIA2+6 *tr*k4 Ft1A2+6*kS5 * tr *tl42-6*)5 *tr+6 *)5 *t1+6 *tr*c45 *tl+... 
6*tr*co6*tl-6*k4 "tl *tr“2-6*k5 “tr-2- 1) /n*2. 

K(5,6)=1/6*(-2*k5 *tr*9 34-2 "kd "3-6" ko *tr*t14243*e56*tr*24+3 *c56*t1°2-67 tr e567 t+... 
Gaks tro 2 *th/he2: 

k(5,10)=- 1/6 *(k4 *tr43-k4 *t1493-3 *k4 *t1*tr424+3 *c45 *tr424+3 *c45 *t1424+-3 *tr*k4 *t142-6*... 
tr*c45*tl)/h‘2; 

k(5,11)=-1/6*(6 *tr*c45 *tl-3 *k5 *tr*t142+3 *k5 *tr*2 *tl+6*tr*c56*tl-k5 *tr43+k5 *t143-... 
3* cS GMO? -3*c56*tI*2-3 *tr*k4 te 3 *k4 At) * tree "13 +k4*te8-3*c45*... 

(O25 ce ot)" 2-675 7tr +675 tno: 

k(5,12)=- 1/6*(k5 *tr*3-k5 *01%3+-3 *c56*tr*24+37c56*tI24-3*k5 *r 102-3 *k5 *tr 2 11-6 *.... 
tr=cso=01)/n*Z: 

KO) = 1/672. kS *tr 34-2 “5 0143-67 k5 *tr “tl 24-3 *c56*tr42-43 e561)" 2-6*tr* C56 tle 
G7K5 the 2 “tle 2; 

k(6,6)=1/6*(2*k5 *tr43-2 *k5 *t143+6*k5 *tr* 142-3 *c56*tr* 2-3 *c56 *tl42-6*)6*tr+6*)6*... 
1+6*tr*c56*tl-6*k5 *tr*2 *tl)/h’2; 

k(6, 1 1)=-1/6*(k5 *tr*3-k5 *tl*3- nBtenGstee24+3 *C56 70124 3 2kS Fel 2-3 *kS tee? *11-... 
6*tr*c56*tl/h’%2; 

K(Ga12 == 17/6 7(-k5 *tr4 3-5 7t143-3 *c56*tr* 2-356 1142-3 *k5 *tr 11 43k) tr 2 Fela 
6*76*tr+6 *j6 *tl+6 *tr*c56*tl)/h’2; 

k@el)=- 1/67 CK 34k) 1433 %C12 *tr42+3*¢12*tl42-3*k] *tr*tl42+3*k 1 tr 2 *tlae 
OF) *tr+67)1*tl-Gel2 str tyne: 

k(7,2)=1/6*(-k1 *tr*3 +k 1 #14343 *c12 *tr*2+3 *c 12*t142-3 *k 1 *tr*t1424+3 *k 1 *tr*2 *tl-... 
6*el2*tr*tbyne 2: 

k(7,7)=1/6*(2*k1 *tr43-2 *k] *143-6*k 1 *tr42*tl+3 *c 12*tr42+3 *c 1 2*1]42-6*) 1 *tr+6*... 
stl 6 kl trl" 2-67 cl 2 ttt) hz: 

k(7,8)=- 1/6" (2 * kl *tr°3-2* 1 "143-6 le tr 2 *tl-3 *c 2 *treerere 12 1142-6 *c12*tr*... 
tl+6*k1 *tr*tl42)/h‘2; 

k(S)1)=1/6*(-k1 *tr*3+k1 *t1434+3*e12*tr424+3*c12*t142-3*K1 *tr*thh24+3*k1 tr 2 *t1-... 
6tel 2 tr thy i*2. 

k(8,2)=- 1/6 *(-3 *k2 *tr*t142+3 *c2 *tr42+3*c2 *t142-6*)2rec(i)*tr+6*)j2rec(i)*tl-6*j2*... 
tr+6*j2*tl-6 *tl *tr*c23-6* tr*c2 *tl+3 *k2*tl*tr*2-6*c 1 2*tr*tl+3 *k1 *tr42*... 
tl-k2*tr43+k2 *t14343 *c23 *tr492+3*c23 *t142-k 1 *tr43 +k 1 *t1434+3*c12 *tr*2+4... 

Seele 142-37 tren. 

Koso) —1/07 (ke? “tn aoe “113-3 TK “ir 12 +3 *c2 3 *tr* 243 *c23 Ftl42+3 *k2 *tl *tr’2-... 
Ol te7e24 /hio2- 

oe i a(n tn 3-2 2k) 3-6 hl i? *t1+3 *c 12 *tr*2+3*cl2 1142-6 *c12*tr™... 
tO 7k 1 *te *tl*2 /h*2* 

k(8,8 )=1/6*(6*k2 *tr*th42+ 3 *c2 *tr42+3 *c2 *tl2-6*)2rec(1)*tr+6*j2rec(1)*tl-6*)j2 *... 


93 


tr+6*j2*tl-6*tl*tr*e23-67 tp 2e2 11-6 P21) *tr"2"6 *cl2 *ir*tl-G*k Ire * 
tl+2*k2 ¥tr*3-2*k2*tl4343*e25tr* 243 C25 117242 7k! *tr*3-2 kK ti 33 
12 *tr424+3*c12*tl*2+6*kl *te*tl* 22: 

k(8,9)=- 1/6*(2*k2 *tr*3-2*k2*t1*3-6*k2 *t)*tr"2+3*e23 *tr*2+3 7023 2-61) 
C25+6°k2 * tr tl 2 ine. 

k(9,2)=1/6*(-k2 *tr43+k2 *t143-3*k2 *tr*tl424+3 *c23 *tr42+3 *c23 *t1424+3 *k2 *t1l *tr42-... 
Gti ste e23)/n 2. 

k(9,3)=- 1/6*(-3 *k2 *tr*t]42-6 *tl *tr*c234+3 *k2 *tl *tr42-3*k3 *tr*t142+3*k3 *tr42 *H1-... 
6*tr*tl *c34-k2 *tr43+k2*tl43+3*c23 *tr424+3*c23*tl42-k3 *tr43+k3 *t143 +... 
3*c34 *tr424+3*c34*1142+3*e3 "tr 243 "C311" 2-6" tl*c3*t-67}3-tit6 45 tl 
6*)3rec(i)*tr+6*j3rec(1)*tl)/h‘?2; 

k(9,4)=1/6 *(-k3 *tr43+k3 *t143-3 *k3 *tr *t142+3 *C34 *tr4243 *034*t1424+3*k3 *tr*2 *tl-... 
6*tr*t]*c34)/h‘2; 

k(9,8)=-1/6*(2 *k2*tr43-2*k2 *t143-6*k2 711 *tr*24+3 *e23 *tr*243*C25 "12-6 ltr. 
c23+6*k2 *tr*tl*2 nee; 

k(9,9)=1/6*(6 *k2 *te *t1*2-6* tl* tr*c23-6*K2 “tl *tr*2+6 *K3 * tr tl” 2-G*k3 “tr 2 *tl- 
6*tr*tl*c34+2 *k2*tr*3-2 *k2 F143 4372 Sete 2432 3 T11% 242 their o- 
2*k3*tl434+3*c34*tr*2+3 *634 *112 43*c3* tr 24s CO tl 2-0 tls tree 
6*}3*tr+6*)3 *tl-6*)3rec(1)*tr+6*j3rec(1)*tl)/h%2; 

k(9, 10)=- 1/6*(2*k3 *tr*3-2*k3 *1148268K3 * tr" t143%eS4* (r°2+ 3 esteem * 
tl*c344+6*k3 *tr*tl*2)/h’2; 

k(10,3)=1/6*(-k3 *tr43+k3*tl43-3*k3 *tr*tl42+3 *c34 *tr424 3 *c34 #14243 *k3 *tr42 ¥ tl-... 
6ttr=t)*c34 )/ee: 

k(10,4)=-1/6*(-6 *tr*c4 *tl-6 *tr*c45 *tl-3 *k3 *tr *t1424+3 *k3 *tr’2 *tl-6* tr *tl *c34-... 
k3*tr43+k3*t1434+3*034 *tr424+3*c34 *t142-3 *tr*k4 *t1424-3*k4 tl *1r42-... 
k4 *tr43.+k4 *t14934-3 *c45 *tr49243*C45 *t1424-3 *c4 *tr42 43 *04 * 1142-674 *trt... 
6*}4*tl-6*)4rec(1)*tr+6*)}4rec(1)*tl)/h%2; 

k(10,5)=1/6*(-k4 *tr43 +k4 *t1493 +3 *k4 *¥ tl ¥tr424+.3*c45 *tr424+.3*c45 *t]42-3 *tr*kK4 *tl42-... 
Ott E45 thin’2: 

k(10,9)=- 1/6*(2*k3 *tr43-2 *k3 *t143-6*k3 *tr*2 *t14+3*c34 *tr424+3*c34 *t142-6* tr*tl*... 
c34+6*k3otr*tl*2 )/noe 

k(10;10)=1/6*(-6*tr*c4 *tl-6*tr*c4'5 *t1+6*k3 *tr*tl“2-6*k3 *tr72*Ul-Gtteedl*c3442~ 
k3 *tr43-2*k3*t1434+3*c34 *tr424+3 #034 *t142 +6 * tr *k4 *t142-6*k4 *tl ¥tr42 4... 
2*k4 *tr*3-2 *k4 * 14343 *c4 5 etr 2 83 * casi 3 *c4 *tr42+3%C4 “e226 " 4 *... 
tr+6*)4 *tl-6 *j4rec(i)*tr+6*)j4rec(1)*tl)/h’2; 

k(10,11)=- 1/6*(2 *k4 *tr*3-2*k4 *t143-6*k4 *1l *0r4%24+3 *c45 *ir* 2294 5 * ee oeomirees4 Sa 
tl+6*tr*k4 *tl42)/h%2; 

k(11,4)=1/6*(-k4 *tr43+k4 *t1434+3*k4 *tl*tr424+3*c45 *tr424+3*c45 142-3 *tr*k4 *t12-... 
6*tr*e454t he, 

k(11,5)=- 1/6*GG6" tr te45 "11-3 *k5 ttr*tl4243*k5 str 2 71-6 tr te5 67 tl-kS tir 3+k5 “4 3+.... 
3*C56*tr42+3 *c5S6*tI42-3 *tr*k4 *t1424+3 *k4 *t] *tr42-k4 *tr493+k4 *t14343 *c45*... 
tr*24+3 *c45 *t142-6*j5 *tr+6*j5 *tl)/h’%2; 

k(11,6)=1/6*k5* tr*34-k 5 "01934-3656 *tr°2+3 *c56711%2-3 "KS *tr th 3 5 tr 2 Ft. 
6*tr*c56*th/A“ 2. 

k(11,10)=-1/6*(2*k4 *tr*3-2 *k4 *t143-6*k4 *tl *tr42+3*c45 *tr424+3*c45 *t142-6*tr*c45 *tl+... 
Gtr k4 Ae ine 

KC1 1,11 )= 164 (2 *k4 *tr*3-2*k4 ¥034-2 Fk5 F0r43-2 *K5 F143 43045 "tr 2435045 T1142 et C56"... 
tr42+3*c56 *tl42-6*k4 *tl*tr42-6*k5 *tr42 *tl-6*j5 *tr+6*)5 *tl+6 *tr*k4 *tl*2+... 
6*k5 *tr*tl42-6* tr*c45 *tl-6 *tr*c56*tl)/h"2; 

k(11,.12)=-1/6*(2*k5 ¥tr*3-2 *k5*t173-6* ko “tr”2 *11+ ike 5Ghtr* 24-3706 112-6 ir C56" 11+. 
G7 ko *tet2y/ ee. 

k(12,5)=1/6*(-k5 *tr*3+k5 *t143+3 *c56 Ftr%2 43 *c5671142-3* ko tr 1192 43 tkS 271... 
Gr i*e5G" lly 2: 

k(12,6)=- 1/6*(-k5 *tr43+k5 *t1434+3*c56*tr424+3 *c56 *t142-3 *k5 *tr*tl424-3 *k5 *tr*2 *tl-... 
6*j6*tr+6 *j6 *tl-6*tr*cS6*tl)/h%2; 
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KOZ, ID== 1/6 *(2 * ko 7143-2 KS 1143-6 She * tre? * 11-3 * 056 * tr 2 +85 6 * 114246 *tr*e56 *tl+... 
Geko st )/h~ 2 

KG 2 = 1/6" (2 koeth 3-2 KS 1h 3-6*k5 *tr* 2 *tl+3*c56 71243 7*c567 t14226*)6*tr+6*j6* .. 
tl+6*k5 *tr*tl42-6*tr*c56*tl )/h’*2; 


Programs used to calculate cylinder gas torques, given angular position data 


at 6, and 6s. 


%o SPRESSr2 

0 - wna n-ne nnn en nnn nnn nnn nnn nn nnn nn nn nn nnn nnn nnn nnn nnn nnn nnn n- 

Jo Program to determine cylinder pressures in 3 cylinder 
Jo two stroke diesel engine, given instantaneous angular 
Jo velocity at two of the six degrees of freedom 

% (crankshaft nose and flywheel) 

Dag ce slap SA ah pe 
ic=[000 -1.2 -1.2 -10.4]*le-04; % set initial conditions 
Gear eee 

% load data 

ae 


load w51k100.md 
t= (0;w5 1k100(1:7938,1)]; 


tt = diff (t); % COMPUTES THE DEL_T’S 

tt = reshape(tt, 126,63): % This phase locks one cycle 

mdt = mean(tt); % Mean dt for each cycle 

[ons co. % Reference each cycle to its own mean 
tt(71) = ttCade/meanitC1)); 

end 


tt = tt* mean(mdt); 

dtrat = tt(1)/mean(tt(2:63,1)); 

teeth = linspace(0,2*p1,127); tooth = 2*pi/127; 

th51 = ic(5) + tooth*dtrat; 

tttt = mean(tt); % This ensemble averages the phases 
cyctm = sum(tttt); % Time for one complete cycle 

thetarS = [ic(5) teeth(1:125)+th51]; 

thetar5( 127) = 2*pitic(5); 

timeS = [0,cumsum(tttt)-(1-dtrat) *tttt(1)]; 

time5(127) = cyctm; 

load wl 1k100.md 

t= (CO wit kl O0Gli7 92051) |: 

t= diff(v); 

tt = reshape(tt,/20,11); % Phase lock one cycle 

t= tei): 

tttt = mean(tt); % Ensemble average the phases 
thetarl = linspace(0,(2*pi),721); % Known positions of the O.E. windows 
time! = [0 cumsum(tttt)]; 

timel = time! *(cyctm/max(time1)); % Adjust times to agree 

omega = (I/cyctm)*2*p1; 

N= ol % set number of nodes for solution 


time = linspace(O,cyctm,N)’; 
dt =cyctm/N; 
theta = linspace(0,2*pi,N)’; 
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fev = 100-dt % filtering cut-off value as ratio of sampling freq 
thetar] = interp|(timel,thetar1,time,’spline); 

thetarS = interp1(time5 ,thetar5,time,’spline); 

thetal = theta + vfilt(thetar1-theta,fcv); 

thetaS = theta + vfilt(thetar5-theta,fcv); 


omegal = deriv(thetal,dt); 
omega5 = deriv(theta5,dt); 
accell = deriv(omegal,dt); 
accel5 = deriv(omega5S,dt); 


figure (1) 

plot (time, deriv(thetar!1,dt), MX’), grid, hold on 

plot (time, omegal, Tr’) 

plot (time, deriv(thetar5 dt), cX’) 

plot (time, omega5, Db’ 

plot (time, ones(size(omega5))*omega, 'g’) 

title (Raw and Filtered Angular Velocity at DOFS | and 5’) 
xlabel (‘time (seconds)’) 

ylabel (angular velocity (rad/sec)’) 

orient landscape 


figure (2) 

plot (time, thetar1!-theta, MX’), grid, hold on 

plot (time, thetal-theta, '’) 

plot (time, thetar5-theta, ‘cX’) 

plot (time, theta5-theta, b’) 

title (Raw and Filtered phase deviation at DOFS | and 5’) 
xlabel (‘time (seconds) 

ylabel (phase deviation from mean (rad)’) 

orient landscape 


1) ne ts oe ee 

% set calibration data 

Ops ctcl eee 

j1=0.02443; Jolb*in*sec*2/rad 
j2=0.2482; J%olb*in*sec*2/rad 
j3=0.1462; %olb*in*sec*2/rad 
j4=0.2482; Jlb*in*sec*2/rad 
[S=1222. %olb*in*sec*2/rad 
j6=0.2870; %olb*in*sec*2/rad 
kKi=3: 1leG: %\|b*in/rad 
k2=7.00e6; %|b*in/rad 
k3=7.00e6; J%\b*in/rad 


k4=10.82e6; %|b*in/rad 
k5=1, 3046; %\b*in/rad 


cl2=00 1 %\b*in*sec/rad 
CZ23=0:0 1 % |b *in*sec/rad 
c34=0.0); % \b*in*sec/rad 
c45=0 01. J \b*in*sec/rad 
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656-001. %|b*in*sec/rad 


c2=0.013; Jolb*in*sec/rad 

Ca=Un ns, %o\b*in*sec/rad 

c4=0.013; J%olb*in*sec/rad 

Tload = 1200; % Ibf*in Load torque 

Fpar = 100; % lof Parasitic force 

Taux = 168; % \bf*in Valvetrain and auxilliary torque 
Tpmp = 9:5; % \bf*in Oil pump torque 

Cetera eee Soe roman Ce 

% Solve equations 6, 5, and 1 

Do BOO ar cia 


global AA BB CC DD TT 
Icomeva= 1 0a.5 | eleinl 1 1]; 
thetar2=zeros(N,1); thetar4=thetar2; thetar6=thetar2; thetar3=thetar2; 
omega2=zeros(N, 1); omegar3=omega2; omega4=omega2; omegar6=omega?2; 
thetar3(1) = 1c(3); 
omegar3(1) = icomega(3); 
Jo------ solve equation 6 for theta6 and omega6------ 
beval = [1c(6),ic(6)+2*p1); 
AA = j0;3BIB =c56; CC = kis 
DD = -Tload + c56*omegaS5 +k5*theta5; 
fif'= zeTrostN), 1); 
kk = zeros(N,N); 
for P= 1-(N-1); 
h = time(i+1)-time(i); 
k = -(AA/h)*[1 -1;-1 1] + (BB/2)*[-1 1;-1 1] + (CC*h/6)*[2 131 2]; 
f = (h/2)*([DDG@);DD(i+1)); 
kkQiit] ji+1) = kkQiit+1,1:1+1) + k; 
ff(iit+l) = ff(isi+l1) + f; 
end 
Jo----- apply boundary conditions ----- 
kk(1,:) = zeros(1,N); Kk(N,:) = zeros(1,N); 
kk(1,1)=1; kk(N,N) = 1; 
ffl = beval(l). 11) =veval(2); 
Jo----- solve matrix eqn ----- 
theta6 = kk\ff; 
omega6 = deriv(theta6,dt); 
Jo------ solve equation 5 for theta4 ------ 
thetar4 = (Taux + k5*(theta5-theta6) + c56*(omega5-omega6) +... 
k4*thetaS + j5*accel5)./k4; 
theta4 = theta + vfilt(thetar4-theta,fcv); 
Jo------ solve equation | for theta2 ------ 
thetar2 = (Tpmp + k1 *thetar] + j1*accell)./k1; 
theta2 = theta + vfilt(thetar2-theta,fcv); 
Yo------ compute omega and accel for dofs 2 and 4 ------ 
omega2 = deriv(theta2,dt); % raw velocities 
omega4 = deriv(theta4,dt); 
accel2 = deriv(omega?2,dt); % calculated accelerations 
accel4 = deriv(omega4,dt); 


R25. % Crankshaft eccentricity (in) 
W = 7.556; % reciprocating weight(lbf) 
¢ = 386; % The acceleration of gravity (in/sec’2) 


D7. 


B= 3.605. % Cylinder Bore (in) 

L = 8.80; % Connecting Rod length (in) 
j2rec=(W*R‘%2/(2*g))*(1-cos(2*theta)); Toil in See. 2 
j3rec=(W*R42/(2*g))*(1-cos(2*(theta-4*pi/3))); Jolb*in*sec*2 
j4rec=(W*R42/(2* 2))*(1-cos(2*(theta-2*pi/3))); %lb*in*sec*2 

Tlcyl = zeros(N,1); T2cyl=Tlcyl; T3cyl = Tlcy!; 

pref = 14.706 + 0.00205 *((omega*60/(2*p1))-634); % ref press (psia) 

N1 = min(find(thetal>=(2*pi/3))); % index for TDC Cyl #3 
N2 = min(find(thetal >=(4 *pi/3))); % index for TDC Cyl #2 
accel3 = zeros(size(accel2)); 

sl] = R*¥cos(theta) + sqrt(L%*2 - (R*2)*sin(theta).%2); 

s2 = R*cos(theta-4*pi/3) + sqrt(L%2 - (R*2)*sin(theta-4*pi/3).%2); 

s3 = R*cos(theta-2*pi/3) + sqrt(L’2 - (R“2)*sin(theta-2*pi/3).%2); 


Sp1 = deriv(sl,dt); Spdl = deriv(Sp1,dt); % Piston speed (in/sec) and 
Sp2 = deriv(s2,dt); Spd2 = deriv(Sp2,dt); % Piston acceleration (in/sec’2) 
Sp3 = deriv(s3,dt); Spd3 = deriv(Sp3,dt); 

Tlrec = -(W/g)*Spd1.*(Sp1/omega); % in*lbf Reciprocating torque 


T2rec = -(W/g)*Spd2.*(Sp2/omega); 

T3rec = -(W/g)*Spd3.*(Sp3/omega); 

Tparl = Fpar*abs(Sp1/omega); % in*\lbf Parasitic torque 
Tpar2 = Fpar*abs(Sp2/omega); 

Tpar3 = Fpar*abs(Sp3/omega); 


Jo------ Step one: Determine known values of Tcyl from pref ------ 
Jo------ (known values of Tcyl are 0) ------ 

Jo------ Step two: Solve for theta3 throughout cycle ------ 

Jo------ solve equation 2 for theta3 ------ 


thetar3(N1:(N2-1)) = (G2+)2rec(N1:(N2-1))).*accel2(N1:(N2-1)) + ... 
c12*(omega2(N1:(N2-1))-omegal(N1:(N2-1))) + ... 
k1*(theta2(N 1:(N2-1))-thetal(N1:(N2-1))) + k2*theta2(N1:(N2-1)) +... 
c2*omega2(N 1:(N2-1)) - Tlcyl(N1:(N2-1)) - Tlrec(N1:(N2-1)) +... 
Tparl(N1:(N2-1)))./k2; 
Jo------ solve equation 3 for theta3 and omega3 ------ 
BB = c234+c344c3; CC = k2+k3; 
DDT = T2cyl + T2rec - Tpar2 + c23*omega2 + k2*theta2 + c34*omega4 + k3*theta4; 
for 1 = 1:(N1-2); 
AA = j3 + j3rec(i:i+1); TT = time(isit+1); DD = DDT(isi4+1); 
[T,X] = ode45(’seqns2’,[time(1),time(i+ 1 )],[thetar3(1);omegar3(1)]); 
thetar3(i+1) = X(length(T), 1); 
omegar3(i+1) = X(length(T),2); 
end 
Jo------ solve equation 4 for theta3 ------ 
thetar3(N2:N) = (G44+j4rec(N2:N)).*accel4(N2:N) + k3*theta4(N2:N) + ... 
(c45+c4)*omega4(N2:N) - c45*omega5(N2:N) + k4*(theta4(N2:N) - ... 
theta5(N2:N)) - T3cyl(N2:N) - T3rec(N2:N) + Tpar3(N2:N))./k3; 
Yo------ filter theta3 and derive omega3 and accel3 ------ 
theta3 = theta + vfilt(thetar3-theta,fcv); 
omega3 = deriv(theta3,dt); 
accel3 = deriv(omega3,dt); 
Jo------ Step three: solve for remaining Tcyl values ------ 
Jo------ solve equation 2 for TIcyl ------ 


Ticyl(1:(N1-1)) =-Tlrec(1:(N1-1)) + Tparl(1:(N1-1)) + G24j2rec(1:(N1-1))).¥accel2(1:(N1-1)) + ... 


c12*(omega2(1:(N1-1))-omegal(1:(N1-1))) + k1*(theta2(1:(N1-1))-thetal (1:(N1-1))) +... 
c23*(omega2(1:(N1-1))-omega3(1:(N1-1))) + k2*(theta2(1:(N1-1))-theta3(1:(N1-1))) +... 
c2*omega2(1:(N1-1)); 

T1cyl(N2:N) = -Tlrec(N2:N) + Tparl(N2:N) + (j24+j2rec(N2:N)).*accel2(N2:N) + ... 
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c12*(omega2(N2:N)-omegal(N2:N)) + k1 *(theta2(N2:N)-thetal(N2:N)) + ... 
c23*(omega2(N2:N)-omega3(N2:N)) + k2*(theta2(N2:N)-theta3(N2:N)) + ... 
c2*omega2(N2:N); 

Jo------ solve equation 3 for T2cyl ------ 

T2cyl(N1:N) = -T2rec(N1:N) + Tpar2(N1:N) + G3+j3rec(N1:N)).*accel3(N1:N) + ... 
c23*(omega3(N1:N)-omega2(N1:N)) + k2*(theta3(N1:N)-theta2(N1:N)) + ... 
c34*(omega3(N1:N)-omega4(N1:N)) + k3*(theta3(N1:N)-theta4(N1:N)) + ... 
c3*omega3(N1:N); 

Jo------ solve equation 4 for T3cyl ------ 


T3cyl(1:(N2-1)) = -T3rec(1:(N2-1)) + Tpar3(1:(N2-1)) + G4+j4rec(1:(N2-1))).*accel4(1:(N2-1)) + ... 


c34*(omega4(1:(N2-1))-omega3(1:(N2-1))) + k3*(theta4(1:(N2-1))-theta3(1:(N2-1))) +... 
c45*(omega4(1:(N2-1))-omega5(1:(N2-1))) + k4*(theta4(1:(N2-1))-theta5(1:(N2-1))) +... 
c4*omega4(1:(N2-1)); 

Jo------ Convert to FI *LBF ------ 

Tiheyl=T lewis: 

T2cyl = T2cyl./12; 

T3cyl = T3cyl./12; 


a 5 (rc 
J Determine and plot measured/predicted cylinder torques 
Gf cd Sn ee aa SES ee OEE os co 
load walk100.md % ECA cylinder #1 pressure data 
load wb1k100.md % ECA cylinder #2 pressure data 
load wclk100.md Jo ECA cylinder #3 pressure data 


pa = reshape (walk100,5,720); Plcyl = [pa(1,:) pa(1)]; 

pb = reshape (wb1k100,5,720); P2cyl = [pb(1,:) pb(1)]; 

pc = reshape (wc1k100,5,720); P3cyl = {pc(1,:) pce(1)]; 

pos = linspace(0,2*pi,721); 

dtm = (2*pi/omega)/720; 

slm = R*cos(pos) + sqrt(L*2 - (R*2)*sin(pos).*2); 

s2m = R*cos(pos-4*pi/3) + sqrt(L*2 - (R*2)*sin(pos-4*pi/3).%2); 
s3m = R*cos(pos-2*pi/3) + sqrt(L’*2 - (R*2)*sin(pos-2*pi/3).%2); 
Splm = deriv(slm,dtm); % Piston speed (in/sec) 

Sp2m = deriv(s2m,dtm); 

Sp3m = deriv(s3m,dtm); 

Tlcylm = -((Plcyl.*pref-pref)*(pi*B%2/4).*Sp1m/omega)/12; % (ft*lbf) 
T2cylm = -((P2cyl.*pref-pref)*(pi*B%*2/4).*Sp2m/omega)/12; % (ft*lbf) 
T3cylm = -((P3cyl.*pref-pref)*(pi*B%2/4).*Sp3m/omega)/12; % (ft*lbf) 
figure(3) 

plot(pos,T1lcylm, bX’,pos,T2cylm, bO’.pos,T3cylm, b.) 

legend(‘Cyl #1’, ‘Cyl #2’,'Cy! #37) 

grid, hold on 

plot(thetal , Tl cyl, rX’,thetal ,T2cyl, rO’»,thetal .T3cyl,r.’) 

title( Measured (blue) and Predicted (red) Cylinder Gas Torques’) 
ylabel( Torque (ft*lbf) 

xlabel(Crank Angle’) 

orient landscape 

figure (4) 

plot(pos,T1lcylm+T2cylm+T3cylm, b’) 

erid, hold on 

plot(thetal ,Tlcyl+T2cyl+T3cyl,r) 

title( Measured (blue) and Predicted (red) Total Cylinder Gas Torque’? 
ylabel (Torque (ft*lbf)) 

xlabel (Crank Angle’) 

orient landscape 

figure(5) 


2) 


degm = 180*pos/p1; 

degp = 180*theta/pi; 

subplot(3,1,1) 

plot(degm,T1cylm,k.’) 

grid, hold on 

plot(degp,TI1cyl,k’) 

legend (Meas ‘Pred 

ylabel (Cyl #1 Torque (ft*lbf)’) 
title’1000 RPM, 100 ft*lbf Cylinder Gas Torques’) 
axis({0,360,-500,1000}) 

subplot(3,1,2) 

plot(degm,T2cylm, k.’) 

grid, hold on 

plot(degp,T2cyl,k’) 

ylabel (Cyl #2 Torque (ft*lbf)’) 
axis({0,360,-500, 1000]) 

subplot(3,1,3) 

plot(degm, T3cylm, k.) 

grid, hold on 

plot(degp,T3cyl,k’) 

ylabel (Cyl#3 Torque (ft*lbf)’) 
xlabel(’‘Crank Angle (deg)’) 
axis((0,360,-500,1000]) 

orient tall 

figure (6) 
plot(degm,T1cylm+T2cylm+T3cylm, k.’) 
grid, hold on 

plot(degp, Tlcyl+T2cyl+T3cyl,k ) 

title (1000 RPM, 100 ft*lbf Total Cylinder Gas Torque’) 
legend(Meas ’°Pred 

ylabel (Total Cylinder Gas Torque (ft*Ibf)?) 
xlabel (Crank Angle (deg)) 
axis((0,360,-400,800]) 


orient tall 
% DERIV 
% Function to determine |-D derivative of a vector using 
Jo a central difference technique. 
% Xd = Deriv(X,dt) 
Jo Returns the derivative of the vector X as a function of 
Jo t, given the time step dt. Default value for dt is 1. 
function xd = deriv(x,dt) 
if nargin == 1; 
d= ae 
end 


N = length(x); 

xd = zeros(size(x)); 

xdf = diti(x); 

xd(2:(N-1)) = (xdf(2:(N-1)) + xdf(1:(N-2)))./(2* dt); 
xd(1) = xdidl dt: 

xd(N) = xdf(N-1 )/dt; 

% Forward difference format 

%xd(1:N-1) = xdf/dt; 

Joxd(N )=xd(1); 
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To SEQNS2 


Jo Function to solve arbitrary 2nd order ode in form 
Jo AAxdd + BBxd + CCx = DD 
Jo where AA, BB, CC, and DD are global variables 


function xdot=seqns2(t,x) 

global AA BB CC DD TT 

DDS = DD(1) + diff(DD)*(t-TT(1))/diff(TT); 
AAS = AA(1) + diff(AA)*(t-TT(1))/diff(TT); 
KOOt =iZerOs( 241). 

XCOUCI ae 

xdot(2) = (DDS - BB.*x(2) - CC.*x(1))./AAS; 


Jo VEE 

Jo Function performs fast fourier transform filtration 

%o of high frequency components of given data 

Jo ia— VEC rey) 

%o Y is filtered data 

Jo X 1s input data 

Jo FCV is frequency cutoff value (as a fraction of 

Jo the sampling frequency) 

function Y = vfilt(x,fcv) 

N = length(x); % number of data points 
re Tih): % Fourier transform of data 
ico = fixgicv*N): % index in D 

DF = zeros(1,N); 

DF(1:(ico+1)) = D(1:(ico+1)); % filter high frequencies 
DF((N-ico+1):N) = D((N-ico+1):N); % mirror image data 

te ete). % inverse fft 

Y = real(Y); 


¥ — reshape(vssize(x)); 


Program used for analysis of torsional model natural frequencies and modal 


shapes 

Of aoe Nee een 

Jo NATFREQ 

Ore ine co Soames eo eee 

Jo Determines natural frequencies of torsional model 
j1=0.02443; Jolb*in*sec*2/rad 
j2=0.2482; Jolb*in*sec*2/rad 
j3=0.1462; Jolb*in*sec*2/rad 
j4=0.2482; Jolb*in*sec*2/rad 
1e=)2 2 Jolb*in*sec*2/rad 
j6=0.2870; %\b*in*sec*2/rad 
jrec=0.02478; 

k1=3.11e6; Jolb*in/rad 
k2=7.00e6; %|b*in/rad 
ko—7/.0Weo, Jo\b*in/rad 


k4=10.82e6; %\b*in/rad 
k5=304e6: %\b*in/rad 


el 2=0 O01 % \b*in*sec/rad 
c23=0:015 %\b*in*sec/rad 
c34=0.01; %\|b*in*sec/rad 
€45=0:01; %|b*in*sec/rad 


10] 


c56=0:01- J%olb*in*sec/rad 


e2=0:0113: %\b*in*sec/rad 
63=0:013: %olb*in*sec/rad 
e4—0:013- J%lb*in*sec/rad 
eo NE es ieee ree 2 re 


J=[j100000;0 j2+jrec 0 0 0 0;0 0 j3+jrec 0 0 0;... 
000 j4tjrec 00;0000 j5 0;00 00 0 j6]; 


K = [kl] -k1 000 0;-k1 k1+k2 -k2 0 0 0;0 -k2 k2+k3 -k3 0 0;... 
O 0 -k3 k3+k4 -k4 0;0 0 O -k4 k4+k5 -k5;0 0 0 O -k5 kS]; 


C7, ee eee pe P nent POSS oo oa. bc ccndeuee eee eee eee 
ws = eig(K/J); 
W = sqrt(ws) % (rad/sec) natural frequencies 
whz = w/(2*pi) % (Hz) natural frequencies 
for 1= 1:6: 
kj = K-J*w(i)’2; 
ACh) = 1; 


A(2,i) = -A(1,1)*kj(1, I)/kjy 1,2); 
A(3,1) = (-A(1,1)*kj(2,1)-A(2,1)*kj(2,2))/ky (2,3); 
A(4,1) = (-A(2,1)*kj(3,2)-A(3,1)* kj (3,3) WKIG,4); 
A(5,1) = (-A(3,1)*kj(4,3)-A(4,1)*kj(4,4))/kj(4,5); 
A(6,1) = (-A(4,1)*kj(5,4)-A(5,1)* kj (5,5) )/kj(5,6); 
end 
[whz,I] = sort(abs(whz)); 
A 
figure (1) 
for 1 = 276; 
subplot(S, 1,1-1) 
plot(A(:,IG)),k), grid, hold on 
plot(A(:,IG)), ko) 
ylabel ({num2str(whz(1),4),’ Hz) 
%Jtitle (( Mode ’.num2str(i, 1)]) 
end 
orient tall 
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APPENDIX F. UNCERTAINTY ANALYSIS 


Various sources of error will affect the accuracy of results in this study. The 
following discussion contains a qualitative summary of the potential sources of 
experimental uncertainty. 

i* Optical Encoder 

The use of the flexible coupling results in a torsional oscillation of the encoder 
disk about the actual angular position of the crankshaft nose (6;). The frequency of the 
oscillation is apparent from the measured data and discussed 1n Appendix C. Heidenhain 
[Ref 21] lists a kinematic error of transfer of +10’ for the encoder coupling, 
corresponding to a vibrational amplitude of about 1 x 10” radians. The oscillation seen 
in experimental data varies from about that value up to 1 x 10° radians. But this high 
frequency oscillation is easily filtered out of the raw data, and its amplitude is still an 
order of magnitude smaller than the amplitude of the crankshaft angular velocity — 
fluctuations. 

The 3600 count optical encoder allows a theoretical adjustment to within 0.1°, but 
this 1s only as accurate as the TDC alignment for the engine. Due to an inadequate 
method of determining TDC, this error may grow to 1 or 2°. 

Te Flywheel 

The TDC position for the flywheel data is determined by a comparison of the first 
At to subsequent values, assuming that the crankshaft has zero twist at TDC. This is a 
reasonable assumption, but not completely accurate. Therefore, an expected error 1s 


introduced due to ambiguity in determination of the flywheel TDC angular position. 
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Tooth to tooth variation will result in a cyclic error for the flywheel data. A 
motoring measurement of the flywheel at a constant speed could be used to generate a 
correction signal, eliminating this error [Ref 8]. This was not done for this study, 
assuming the error would be smal] relative to the measured velocity fluctuations. 

A radial vibration is present in the flywheel during engine operation. Since the 
proximeter is mounted to view the ring gear teeth from the side, this radial vibration will 
cause the position of the proximeter relative to the tooth to oscillate up and down. This 
up and down oscillation will induce a cyclic variation in the pulse width not associated 
with crankshaft rotational velocity fluctuations. However, it is expected that this 
Variation is reasonably small enough to be ignored. 

> Measurement Error 

Engine speed and load vary slowly during data collection. Since data for the 
various runs are not collected simultaneously, there are small differences 1n speed and 
load for data comprising a set. Typically, engine speed variation was within 20 RPM of 
the stated value, and load was within | ft*lbf. A correction is made in the pressure 


prediction program to account for small differences in rotation speed between the 9, and 


Qs; data. 
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